spectractor.serve

1 2 '''Spectractor submodule containing functions for managing various files, 3 journals, creating logs etc. 4 5 @author: Dmitry Nasonov 6 ''' 7 8 import re 9 import os, sys, shutil, glob 10 import cPickle 11 import collections 12 13 from spectractor import get_wv_range 14 data_dir, out_dir = 'data', 'out' 15 16

17 -def filler(indct, key, val):

18 """Fill dictionary which values are lists.""" 19 try: 20 if val not in indct[key]: 21 indct[key].append(val) 22 except KeyError: 23 indct[key] = [val, ] 24 except TypeError: 25 print >> sys.stderr, "TypeError in filler! Key, value:", key, val 26 sys.exit(True)

27 28

29 -def fill_dctset(indct, key, val):

30 """Fill dictionary which values are sets.""" 31 try: 32 indct[key].add(val) 33 except KeyError: 34 indct[key] = set([val])

35 36

37 -def backuper(file_pth, verbose=False):

38 """Backup file if it is exist: move it with '-backup' postfix.""" 39 if os.path.isfile(file_pth): 40 shutil.move(file_pth, file_pth + '-backup') 41 if verbose: 42 print "Move", file_pth, "to", file_pth + '-backup'

43 44

45 -def rexper(regexp, path, reflag=None, splitnig=False):

46 """Get parameters from file name using regexp. 47 Parameters are: int ID (for example, night number), flag. 48 For now, regexp.findall must return 3 values! If not, return (None, None). 49 @param regexp: regular expression describing file name 50 @param path: path to file 51 @param reflag: tuple of two values. If not empty, flag containing in first 52 value will be replaced to second one 53 @param splitnig: split ID to tuple of 'night' and spectrum number 54 @return: ID (night number): int or tuple of ints with splitnig, flag 55 """ 56 try: 57 bname = os.path.basename(path) 58 night, img, flag = regexp.findall(bname)[0] 59 except (IndexError, ValueError): 60 print >> sys.stderr, "Filename", bname, "does not match regexp!" 61 return None, None 62 if reflag and flag in reflag[0]: 63 flag = reflag[1] 64 if not splitnig: 65 night = int(night + img.zfill(3)) 66 else: 67 night = (int(night), int(img)) 68 return night, flag

69 70

71 -class Walker:

72 """Directory walker: search files and dirs in args matching regexp. 73 @param exts: set of file extensions for matching 74 @param namrxp: filename regular expression, but without extension 75 @param args: list of paths to files and dirs 76 """

77 - def __init__(self, args, exts=("200", "fds"), namrxp=".*", verbose=False):

78 self.files_dct = {} 79 self.workdirs = [] 80 self.path_set, self.spath_set = set(), set() 81 self.verbose = verbose 82 # rexp_spec contains 3 groups: night number, spectrum number and flag. 83 self.rexp_spec = "^[a-z_]{0,3}(\d{3})(\d{1,3})[-_]?([a-z_\-]*).*\." 84 self.fdsrexp = re.compile(self.rexp_spec+'fds'+'$', re.I | re.U) 85 # rexp_mt and rexp_fits contain groups for date (year, month, day), 86 # spectrum number and flag. 87 self.rexp_mt = \ 88 "^[lns](\d{1})([A-C0-9])(\d{2})[_-]?(\d{2,3})([a-z_]?).*" 89 self.rexp_fits = \ 90 "^(?:Bn)?(\d{4})(\d{2})(\d{2})[_-]?(\d{3})([spbfdot_]?).*" 91 exts = "|".join(exts) 92 self.rxp = re.compile(namrxp+'\.('+exts+')$', re.I | re.U) 93 # Pathfinder: find dirs and files in args, append them. 94 for arg in args: 95 if os.path.isdir(arg): 96 self.workdirs.append(arg) 97 if verbose: 98 print "Appending", arg, "to the list of working dirs" 99 elif os.path.isdir(os.path.join(os.getcwd(), arg)): 100 self.workdirs.append(os.path.join(os.getcwd(), arg)) 101 if verbose: 102 print "Appending", os.path.join(os.getcwd(), arg), \ 103 "to the list of working dirs" 104 elif os.path.isfile(arg) and len(self.rxp.findall(arg)): 105 filler(self.files_dct, self.rxp.findall(arg)[0], arg) 106 if verbose: 107 print "Appending", arg, "to the list of working files" 108 elif verbose: 109 print arg, "is not a directory or needed file, continue."

110

111 - def dirwalker(self):

112 """Walk on workdirs, collect subdirectories. 113 @return: list of subdirectories 114 """ 115 if not self.workdirs: 116 print "Sorry, no dirs for walk." 117 return 118 subdirs = [] 119 for dirw in self.workdirs: 120 # dir path, inner dirs, inner filenames 121 for root, dirs, files in os.walk(dirw): 122 for wdir in sorted(dirs): 123 subdirs.append(os.path.join(root, wdir)) 124 return subdirs

125

126 - def walk(self, exclude="excluded"):

127 """Walk into working dirs, search some kind of files. 128 @param exclude: if it in file path, exclude file. Useful when 129 we want to exclude some subdir 130 """ 131 if not self.workdirs: 132 if not self.files_dct: 133 print "Sorry, no files found." 134 return 135 for dirw in self.workdirs: 136 # dir path, inner dirs, inner filenames 137 for root, dirs, files in os.walk(dirw): 138 #sorted(dirs + files) == sorted(os.listdir(root)) 139 for wfile in sorted(files): 140 if len(self.rxp.findall(wfile)) and exclude not in root: 141 filler(self.files_dct, self.rxp.findall(wfile)[0], 142 os.path.join(root, wfile))

143

144 - def resolver(self, ext='rv', filtr=None, reflag=None, splitnig=False):

145 """Select data, get ID and 'flag' from filename, collect paths. 146 Here we assume, that file name contains information about ID 147 (night number, number of spectrum) and flag. Others will be broken 148 (or you must change default regular expression). 149 @param ext: select file type, must be file extension 150 @param filtr: selecting only given flags, exclude others. Select all 151 when filtr in ("all", None) 152 @param reflag: tuple of two values. If not empty, flag containing in 153 first value will be replaced to second one 154 @param splitnig: split ID to tuple of 'night' and spectrum number 155 @return: list of tuples with 'night', flag and path to file 156 """ 157 try: 158 files = self.files_dct[ext] 159 except KeyError: 160 print >> sys.stderr, "Sorry,", ext, "files seem not to exist." 161 return self.path_set 162 rexp = re.compile(self.rexp_spec + ext + '$', re.I | re.U) 163 for path in files: 164 night, flag = rexper(rexp, path, reflag=reflag, splitnig=splitnig) 165 # Collect "kosher" paths: 166 if filtr in ("all", None) or flag in filtr: 167 self.path_set.add((night, flag, path)) 168 elif self.verbose: 169 print "Exclude", night, path 170 if not len(self.path_set): 171 print >> sys.stderr, "No", ext, "files in path!" 172 return self.path_set

173

174 - def select_spectra(self, journal, ext='fits', iflag=None, splitnig=False):

175 """Select data, get ID ('night') and flag from filename, collect paths. 176 Here we assume, that file name contains information about ID 177 (night number, number of spectrum) and flag. Others will be broken 178 (or you must change default regular expression). Then we check presence 179 of ID in journal keys, get spectrum parameters from journal and append 180 this data to resulting list. 181 @param journal: journal dictionary with spectra IDs as keys and spectra 182 parameters (e.g. JYear, MJD, spectrograph key, limiting 183 wavelengths, S/N etc) as values 184 @param ext: select file type, must be file extension 185 @param iflag: Flags to select. Can be a string, if flag is one symbol 186 @param splitnig: split 'night' to tuple of night and spectrum number 187 @return: list of tuples with ID ('night'), flag, path to file and 188 other parameters. Currently parameters are fds path, dictionary 189 of points, defining the flatten polynome, heliocentric 190 correction and S/N value 191 """ 192 try: 193 files = self.files_dct[ext] 194 except KeyError: 195 print >> sys.stderr, "Sorry,", ext, "files seem not to exist" 196 return self.spath_set 197 rexp = re.compile(self.rexp_spec + ext + '$', re.I | re.U) 198 if not journal: 199 self.journal = self.mk_dummy_journal(rexp, files) 200 else: 201 self.journal = journal 202 Sp = collections.namedtuple("Sp", "id flag pth fdpth plpth Vr sn jy jd") 203 for path in files: 204 night, flag = rexper(rexp, path, reflag=False, splitnig=splitnig) 205 if flag and iflag and flag not in iflag: 206 if self.verbose: 207 print "Exclude", night, "with flag", flag 208 continue 209 if night not in self.journal: 210 continue 211 jy, mjd, spc, lwv, hwv, Vr, sn, fdsnam = self.journal[night][:8] 212 #jy, mjd, spc, lwv, hwv, Vr, sn, fdsnum = self.journal[night] 213 pls_pth = os.path.splitext(path)[0]+'.ccm.txt' 214 if not os.path.isfile(pls_pth): 215 print >> sys.stderr, "Exclude", path+": no ccm.txt file" 216 continue 217 fds_pth = os.path.join(os.path.dirname(path), str(fdsnam)+'.fds') 218 if not os.path.isfile(fds_pth): 219 fds_pth = self.find_fds(path, str(night).zfill(6)) 220 if not fds_pth: 221 print >> sys.stderr, "Exclude", path+": no fds file" 222 continue 223 if self.verbose: 224 print "Append", night, flag, "with shift", Vr 225 params = Sp(night, flag, path, fds_pth, pls_pth, Vr, sn, jy, mjd) 226 self.spath_set.add(params) 227 return self.spath_set

228

229 - def mk_dummy_journal(self, rexp, files):

230 """Make dummy journal with paths to fds and wavelength limits 231 """ 232 journal = {} 233 for path in files: 234 night, flag = rexper(rexp, path, reflag=False, splitnig=False) 235 if night not in journal: 236 fds_pth = self.find_fds(path, str(night).zfill(6)) 237 if not fds_pth: 238 print "Exclude", path, "no fds" 239 continue 240 lwv, hwv = get_wv_range(fds_pth) 241 fds_pth = os.path.basename(fds_pth) 242 #jy, mjd, spc, lwv, hwv, Vr, sn, fdsnam 243 journal[night] = [0, 0, '', lwv, hwv, 0, 1, fds_pth] 244 return journal

245

246 - def find_fds(self, path, night):

247 """Try to find fds, using path to file with spectrum. 248 We assume that fds and spectrum name are similar and search fds 249 given at the same night, that of spectrum. 250 """ 251 fdses = {} 252 nig, num = night[:3], int(night[3:]) 253 rxp = '*'+nig+'*'+'.fds' 254 if os.path.dirname(path): 255 rxp = os.path.dirname(path) + os.sep + rxp 256 fds_pths = glob.glob(rxp) 257 if not fds_pths: 258 return 259 for fdpth in fds_pths: 260 fdnam = os.path.basename(fdpth) 261 night, img, flag = self.fdsrexp.findall(fdnam)[0] 262 if img: 263 fdses[abs(num-int(img))] = fdpth 264 return fdses[min(fdses)]

265

266 - def select_images(self):

267 """Select image files with filenames matching rexp_fits or rexp_mt. 268 Currently only fits, fts and mt file extensions are supported. 269 """ 270 rexpfts = re.compile(self.rexp_fits, re.I | re.U) 271 rexpmt = re.compile(self.rexp_mt, re.I | re.U) 272 Sp = collections.namedtuple("Sp", "pth num date flag keymode") 273 for ext in self.files_dct: 274 for pth in self.files_dct[ext]: 275 fnam = os.path.basename(pth) 276 if ext.lower() in ("fits", "fts") \ 277 and len(rexpfts.findall(fnam)): 278 year, mon, day, num, flag = rexpfts.findall(fnam)[0] 279 elif ext.lower() == "mt" and len(rexpmt.findall(fnam)): 280 year, mon, day, num, flag = rexpmt.findall(fnam)[0] 281 mon = str(int(mon, 16)).zfill(2) 282 year = str(2000 + int(year)) if int(year) else "2010" 283 else: 284 continue 285 keymode = False 286 # Set flag for incorrect FITS header cards in NES spectra 287 # taken in 2008: 288 if 20080422 < int("".join((year, mon, day))) < 20081103: 289 keymode = "flstr" 290 date = tuple(map(int, (year, mon, day))) 291 params = Sp(pth, num, date, flag, keymode) 292 self.path_set.add(params) 293 return self.path_set

294 295

296 -def read_journal(j_pth): return cPickle.load(open(j_pth, 'rb'))

297 298

299 -def mod_journal(journal, sn=10, rvc=True, rvcorr_dct=None, instr=None, 300 lams=None, jd=None, verbose=False):

301 """Exclude or correct items in journal depending on input parameters. 302 @param journal: dictionary with spectra main parameters 303 @param sn: limiting signal-to-noise ratio value: exclude all spectra with 304 lower S/N 305 @param rvc: Set Va shift to zero. Useful for telluric line measurements 306 @param rvcorr_dct: Dictionary of heliocentric velocity corrections, 307 defined, for example, by telluric lines measurements 308 @param instr: Select spectra obtained with given spectrograph 309 @param lams: Exclude all spectra, which no contain given vawelength 310 @param jd: Exclude all spectra, which obtained earlier, than given JYear 311 value 312 """ 313 if not rvc: 314 print "Set Va to zero" 315 for nig, dat in journal.items(): 316 jy, mjd, spc, lwv, hwv, Vr, jsn, fn = dat[:8] 317 if sn and jsn < sn: 318 if verbose: 319 print "Exclude", nig, "S/N", sn 320 del journal[nig] 321 continue 322 if instr and spc not in instr: 323 if verbose: 324 print "Exclude", nig, "instr", instr 325 del journal[nig] 326 continue 327 if lams and (lwv > lams[1] or hwv < lams[0]): 328 if verbose: 329 print "Exclude", nig, "wv limits", lwv, hwv 330 del journal[nig] 331 continue 332 if jd and jy < jd: 333 if verbose: 334 print "Exclude", nig, "JD", jd 335 del journal[nig] 336 continue 337 if not rvc: 338 if type(journal[nig]) is not list: 339 journal[nig] = list(journal[nig]) 340 journal[nig][5] = 0 341 if rvcorr_dct and nig in rvcorr_dct: 342 if type(journal[nig]) is not list: 343 journal[nig] = list(journal[nig]) 344 journal[nig][5] += rvcorr_dct[nig]

345 346

347 -def clocker(val, divc=1.):

348 """Convert float value to tuple of (int h, int m, float s). 349 Useful for coordinates and times. 350 @param divc: division coefficient for val 351 """ 352 h, rest = divmod(val/divc, 1) 353 m, rest = divmod(rest * 60, 1) 354 s = round(rest * 60, 2) 355 return int(h), int(m), s

356 357

358 -class Logger:

359 """Extract information from FITS headers of given FITS files, create log 360 @param fitsfixmode: when FIX header is incorrect, try to fix it with given 361 option. Can be 'fix' and 'silentfix' 362 """

363 - def __init__(self, fitsfixmode='silentfix', verbose=False):

364 self.fitsfixmode = fitsfixmode 365 self.verbose = verbose 366 # Define keys, presented in FITS headers: 367 self.stdkeys = ['AUTHOR', 'REFERENC', 'DATE', 'ORIGIN', 'DATE-OBS', 368 'TELESCOP', 'INSTRUME', 'OBSERVER', 'OBJECT', 'EQUINOX', 'EPOCH'] 369 self.nonstdkeys = ['OBSERVAT', 'CHIPID', 'DETNAME', 'CCDTEMP', 370 'P_DEWAR', 'SEEING', 'PROG-ID'] 371 #self.otherkeys = ['WIND', 'T_OUT', 'T_IN', 'T_MIRR', 'PRES', 372 #'ST', 'MT', 'FOCUS'] 373 self.journal_dct = {} 374 self.obs_dct = {} 375 self.key_dct = {"bias": 'b', "thar": 't', "flatfield": 'f', "ff": 'f', 376 "ll": 'f', "sky": 's', "ref.ima.": 'r'} 377 self.postfix, self.ext = '.', 'fits' 378 self.print_str = "[01;31m%s not found in header! Set it to zero.[0m" 379 self.keyset = lambda key: self.key_dct.get(key.lower()) \ 380 if self.key_dct.get(key.lower()) else 'o' 381 self.trygetkey = lambda head, *args: filter(head.has_key, args) 382 self.seper = lambda fstr, h, m, s: \ 383 "%s" % ":".join(["%02d" % h, "%02d" % m, fstr % s]) 384 self.biases, self.flats, self.calib = {}, {}, {}

385

386 - def modkey(self, val, keymode=False):

387 """Modify input value, depending on keymode. 388 Because of different standards of writing FITS headers in SAO RAS, 389 we need to extract the right value. For example, val maybe string, 390 float(string) or int(float * const). 391 @param keymode: False (None), 'flstr' or 'str'. Define, how we can 392 extract correct information from given value 393 """ 394 if val and not keymode and 3000 < abs(val) < 1250000: 395 val = val / 3600. 396 elif keymode == "flstr": # or abs(val) > 1250000: 397 val = str(val) 398 val = float(val[-4:])/3600 + int(val[-6:-4])/60. + int(val[:-6]) 399 elif keymode == "str": 400 if type(val) is str and " " in val: 401 val = val.split() 402 elif type(val) is str and ":" in val: 403 val = val.split(":") 404 val = float(val[2])/3600 + int(val[1])/60. + int(val[0]) 405 return val

406

407 - def pickhead(self, head, keymode, *args):

408 """Look into FITS header for existance of keys, return head[key]. 409 First, we look on existence of keys, containing in args. 410 Second, we get existing value from header or return zero. 411 """ 412 key = self.trygetkey(head, *args) 413 if key: 414 val = head.get(key[0]) 415 else: 416 print self.print_str % " ".join(args) 417 return 0 418 if keymode == "raw": 419 return val 420 else: 421 return self.modkey(val, keymode=keymode)

422

423 - def getcoords(self, head, keymode=False):

424 ra = self.pickhead(head, keymode, 'RA') 425 if not keymode: 426 ra = ra / 15. 427 dec = self.pickhead(head, keymode, 'DEC') 428 az = self.pickhead(head, keymode, 'A_PNT', 'AZIMUTH') 429 zen = self.pickhead(head, keymode, 'Z_PNT', 'ZDIST', 'ZENDIST') 430 if type(az) in (int, float) and az < 0: 431 az += 360 432 #sidertime_start = self.pickhead(head, keymode, 'ST') 433 ra, dec = clocker(ra), clocker(dec) 434 az, zen = clocker(az), clocker(zen) 435 return ra, dec, az, zen

436

437 - def logger(self, curfts, num, date, file_pth, keymode=False):

438 """Create an observation log for input images. 439 Search for some specific keys in given fits header and write it 440 to 'journal': journal_dct and obs_dct dictionaries. 441 """ 442 # Read coords 443 try: 444 ra, dec, az, zen = self.getcoords(curfts[0].header, keymode=keymode) 445 except (ValueError, TypeError, AttributeError): 446 print "Trying to fix", file_pth, "FITS header." 447 curfts.verify(self.fitsfixmode) 448 keymode = "str" 449 ra, dec, az, zen = self.getcoords(curfts[0].header, keymode=keymode) 450 head = curfts[0].header 451 #except (TypeError, AttributeError): 452 #ra, dec, az, zen = (0, 0, 0), (0, 0, 0), (0, 0, 0), (0, 0, 0) 453 # 199xxxxx(DVD01) 20030810(DVD07)--20080331(DVD31), since 20081103_002 454 wind = self.pickhead(head, 'raw', 'WIND') 455 # This block in 20030810(DVD07)--20080331(DVD31), since 20081103_002 456 moscowtime_start = self.pickhead(head, 'raw', 'MT') / 3600. 457 focus = self.pickhead(head, 'raw', 'FOCUS') 458 t_out = self.pickhead(head, 'raw', 'T_OUT') 459 t_in = self.pickhead(head, 'raw', 'T_IN') 460 t_mirr = self.pickhead(head, 'raw', 'T_MIRR') 461 pres = self.pickhead(head, 'raw', 'PRES') 462 # HUMD in 20070902(DVD28)--20080331(DVD31), since 20081103_002 463 hum = self.pickhead(head, 'raw', 'HUMD') 464 # JD since 20080422(DVD31) 465 jd = head.get('JD') 466 # keymode is here! 467 self.journal_dct['KEYMODE'] = [keymode] 468 for key in (self.nonstdkeys + self.stdkeys): 469 if key not in self.journal_dct and head.get(key): 470 self.journal_dct[key] = [head.get(key),] 471 elif head.get(key) and head.get(key) not in self.journal_dct[key]: 472 self.journal_dct[key].append(head.get(key)) 473 elif key not in self.journal_dct and not head.get(key): 474 self.journal_dct[key] = [] 475 # Block of standard "observation" keys 476 time_end = self.pickhead(head, 'raw', 'TM_END') / 3600. 477 try: 478 time_start = self.pickhead(head, 'raw', 'TM-START', 479 'TM_START') / 3600. 480 time = clocker(time_start) 481 except TypeError: 482 time = self.pickhead(head, 'raw', 'TM-START', 'TM_START') 483 exp = int(self.pickhead(head, 'raw', 'EXPTIME')) 484 obj = self.pickhead(head, 'raw', 'OBJECT') 485 flag = self.keyset(obj) 486 if flag == "b" and exp != 0: 487 flag = "o" 488 elif exp == 0 and flag != "b": 489 flag, obj = "b", "BIAS" 490 491 if flag == "b": 492 newname = "".join(('b', num, self.postfix, self.ext)) 493 self.biases[file_pth] = newname 494 elif flag == "f": 495 newname = "".join(('f', num, self.postfix, self.ext)) 496 self.flats[file_pth] = newname 497 elif flag == "t": 498 newname = "".join(('l', num, "t", self.postfix, self.ext)) 499 self.calib[file_pth] = newname 500 501 dateobs = head.get('DATE-OBS') 502 if self.verbose: 503 self.warner(exp, obj, num) 504 self.obs_dct[int(num)] = [date, flag, obj, time, exp, ra, dec, 505 az, zen, dateobs]

506

507 - def warner(self, exp, obj, num):

508 """Warn about some problems.""" 509 if exp == 0 and str(obj).lower() != "bias": 510 print "Wrong object name:", obj, "with exp=0 in spectrum", num

511 #if head.get('CCDTEMP') != None and -103 > head.get('CCDTEMP') > -97: 512 #print "incorrect CCDTEMP", ccdtemp, "in spectrum number", num 513

514 - def logtostr(self):

515 log = [] 516 for num, val in sorted(self.obs_dct.items()): 517 log.append(self.logstringer(num, val)) 518 return log

519

520 - def journaltostr(self):

521 jrn = [] 522 for key, val in sorted(self.journal_dct.items()): 523 if not val or val == [" "]: 524 continue 525 elif key in ("CCDTEMP", "P_DEWAR") and val: 526 val = [str(sum(val)/len(val))] 527 elif key in ("DATE-OBS", "ORIGIN", "AUTHOR", "OBSERVER", 528 "SEEING") and val: 529 val = [", ".join(map(str, val))] 530 jrn.append("%8s %s" % (key, " ".join(map(str, val)))) 531 return jrn

532

533 - def logstringer(self, num, val):

534 """Make a string for given parameters of observation log.""" 535 date, flag, obj, time, exp, ra, dec, az, zen, dateobs = val 536 date = "".join((str(date[0]), str(date[1]).zfill(2), 537 str(date[2]).zfill(2)))+" "+flag 538 #date = "".join(map(str, date))+" "+flag 539 time = self.seper("%02d", *time) 540 ra = self.seper("%04.1f", *ra) 541 dec = self.seper("%04.1f", *dec) 542 az = self.seper("%04.1f", *az) 543 zen = self.seper("%04.1f", *zen) 544 return "%3d %10s %15s %8s %4d %11s %11s %11s %10s %s" % (int(num), 545 date, obj, time, exp, ra, dec, az, zen, dateobs)

546

Source Code for Module spectractor.serve