The VisAO Camera
visaoutils.py
Go to the documentation of this file.
1 ##@file visaoutils.py
2 #
3 # Various utilities for use in VisAOScript
4 #
5 #
6 
7 """@package visaoutils
8 
9  Various utilities for use in VisAOScript
10 
11 """
12 
13 # -*- coding: utf-8 -*-
14 
15 import sys, os, time, calendar, glob, math, copy, pyfits, shutil
16 
17 
18 def get_visao_filename(prefix, t):
19  """
20  Gets the standard VisAO filename (without the leading base name) for a time t.
21  Input:
22  prefix = prefix of the filename
23  t = the time in seconds since the epoch.
24  Returns: a string containing the VisAO filename time portion.
25  """
26  gmt = time.gmtime(t)
27  t_usec = int((t - math.floor(t))*1e6)
28 
29  res = "%s_%04i%02i%02i%02i%02i%02i%06i" % (prefix, gmt.tm_year, gmt.tm_mon, gmt.tm_mday, gmt.tm_hour, gmt.tm_min, gmt.tm_sec, t_usec)
30 
31  return res
32 
33 def get_visao_filename_now(prefix):
34  """
35  Gets a filename
36  """
37  sd = get_visao_filename(prefix, time.time())
38 
39  return sd
40 
41 def visao_wait(waitt):
42  """
43  Signal proof pause
44  """
45  t = time.mktime(time.localtime())
46  while(time.mktime(time.localtime()) - t < waitt):
47  time.sleep(waitt)
48 
49 
50 def sdssFilters():
51  """
52  Return a sequence containing the names of the VisAO bandpasses
53  """
54  filters=["r'", "i'", "z'", "Ys"]
55  return filters
56 
57 def xDither(scale = 1.0):
58  """
59  Return a the x offsets for a standard 5 point dither pattern (starting from 0)
60 
61  scale = multiplicative scale to multiply by, default is 1.0
62  """
63  dx = [-.5*scale, 1.*scale, 0.*scale, -1.*scale]
64 
65  return dx
66 
67 def yDither(scale = 1.0):
68  """
69  Return the y offsets for a standard 5 point dither pattern (starting from 0)
70 
71  scale = multiplicative scale to multiply by, default is 1.0
72  """
73  dy = [.5*scale, 0.*scale, -1.*scale, 0.*scale]
74 
75  return dy
76 
77 def visao_say(txt):
78  """
79  Say something using the MagAO sound server as the VisAO personality.
80 
81  txt = the words to say
82  """
83  cmd = "echo vicki %s | nc zorro.lco.cl 50000" % (txt)
84  os.system(cmd)
85 
86 def visao_script_complete():
87  visao_say("viz ay oh script complete")
88 
89 def visao_script_error():
90  visao_say("viz ay oh script airor")
91 
92 def visao_script_stopped():
93  visao_say("viz ay oh script stopped")
94 
95 def visao_alert(n=1):
96  """
97  Issue an audible alert to the terminal.
98 
99  n = the number of alerts to issue with a 2 second wait between.
100  """
101  if n == 1:
102  print('\a')
103  else:
104  for i in range(n):
105  print('\a')
106  visao_wait(2.)
107 
108 
109 def parseSysLog(fname):
110  f=open(fname, 'r')
111 
112  res = list()
113  for line in f:
114 
115  elem = line.split()
116 
117  tv_sec = float(elem[0])
118  tv_used = float(elem[1])
119 
120  fsec = tv_sec + tv_used/1e6
121 
122  elem[0] = fsec
123  elem.pop(1)
124 
125  res.append(elem)
126  #print elem
127 
128  f.close()
129 
130  return res
131 
132 
133 def timestampFileName(fname, prefix="V47_"):
134  off = len(prefix)
135  yr = int(fname[off:off+4])
136  mo = int(fname[off+4:off+6])
137  d = int(fname[off+6:off+8])
138  hr = int(fname[off+8:off+10])
139  mn = int(fname[off+10:off+12])
140  s = int(fname[off+12:off+14])
141  usec = int(fname[off+14:off+20])
142 
143  tv = calendar.timegm([yr, mo, d, hr, mn, float(s)+float(usec)/1e6, '', '' , ''])
144 
145  return tv
146 
147 
148 def getFileTimes(dir, prefix="V47_", ext=".fits"):
149 
150  if(dir[len(dir)-1] != '/'):
151  dir = dir + '/'
152 
153  fnames = glob.glob(dir+prefix+"*"+ext)
154 
155  res = list()
156 
157  for dfn in fnames:
158  fn = dfn[len(dir):len(dfn)]
159  ft = timestampFileName(fn, prefix)
160  if(ext == ".fits"):
161  joetimes = getJoeTimes(dfn)
162  else:
163  joetimes = [-1,-1]
164 
165  res.append([dfn, fn, ft, joetimes[0], joetimes[1]])
166 
167  return sorted(res, key=lambda times: times[2])
168 
169 def getOneFileTime(dir, prefix="V47_", ext=".fits"):
170 
171  if(dir[len(dir)-1] != '/'):
172  dir = dir + '/'
173 
174  fnames = glob.glob(dir+prefix+"*"+ext)
175  fnames.sort(key=os.path.getmtime)
176 
177  res = list()
178 
179  if len(fnames) < 1:
180  return res
181 
182  dfn = fnames[-1]
183  fn = dfn[len(dir):len(dfn)]
184  ft = timestampFileName(fn, prefix)
185  if(ext == ".fits"):
186  joetimes = getJoeTimes(dfn)
187  else:
188  joetimes = [-1,-1]
189 
190  res.append([dfn, fn, ft, joetimes[0], joetimes[1]])
191 
192  return sorted(res, key=lambda times: times[2])
193 
194 def chooseLogFile(ftime, loglist):
195 
196  for i in range(len(loglist)-1):
197  if(loglist[i][2] <= ftime and loglist[i+1][2] >= ftime):
198  return i
199 
200  return -1
201 
202 def getVisAOExpTimes(ftime, rotime, exptime):
203 
204  expend = ftime-rotime
205  expst = ftime-rotime-exptime
206 
207  return [expst, expend]
208 
209 def getLogEntries(ftime1, ftime2, loglist):
210 
211  #exptime1 = getVisAOExpTimes(ftime1, rotime, exptime)
212 
213  #ist = chooseLogFile(exptime1[0], loglist)
214 
215  ist = chooseLogFile(ftime1, loglist)
216 
217  #exptime2 = getVisAOExpTimes(ftime2, rotime, exptime)
218 
219  #iend = chooseLogFile(exptime2[1], loglist)
220  iend = chooseLogFile(ftime2, loglist)
221 
222  #Just to be safe go one more
223  if(iend < len(loglist) - 1):
224  iend = iend + 1
225 
226  entries = list()
227 
228  entries = parseSysLog(loglist[ist][0])
229 
230  for i in range(ist+1, iend):
231  entries.extend(parseSysLog(loglist[i][0]))
232 
233  return entries
234 
235 
236 def getCCD47MinExp(pxrt,window, bin):
237  """
238  Get the minimum exposure time (also the readout time) for the VisAO CCD47 image.
239  input:
240  pxrt = the pixel rate (2500,250, or 80)
241  window = the window size (1024, 512,256, 64, or 32)
242  bin = the binning (1,2,16)
243  output:
244  minexptime = the minimum exposure time
245  """
246 
247  minexptime = 0
248  if(bin == 1):
249  if(window == 1024):
250  if(pxrt == 2500):
251  minexptime = 1./3.528
252  if(pxrt == 250):
253  minexptime = 1./0.440
254  if(pxrt == 80):
255  minexptime = 1./0.144
256  if(window == 512):
257  if(pxrt == 2500):
258  minexptime = 1./6.703
259  if(pxrt == 250):
260  minexptime = 1./1.487
261  if(pxrt == 80):
262  minexptime = 1./0.535
263  if(window == 256):
264  if(pxrt == 80):
265  minexptime = 1./1.772
266  if(window == 64):
267  if(pxrt == 2500):
268  minexptime = 1./31.283
269  if(window == 32):
270  if(pxrt == 2500):
271  minexptime = 1./42.779
272  if(bin == 2):
273  if(window==1024):
274  if(pxrt==80):
275  minexptime = 1./0.551
276  if(bin==16):
277  if(window==1024):
278  if(pxrt==80):
279  minexptime=1./10.420
280 
281  if(minexptime == 0):
282  print 'getCCD47MinExp: Read out time not found for parameters'
283  return -1
284 
285  return minexptime
286 
287 
288 def getJoeTimes(fname):
289  """
290  Get the exposure time and readout time for a VisAO CCD47 image.
291  input:
292  fname = file name
293  output:
294  [rotime, exptime]
295  """
296  hdu = pyfits.open(fname)
297 
298  exptime = float(hdu[0].header["EXPTIME"])
299  pxrt = float(hdu[0].header["V47PIXRT"])
300  window = float(hdu[0].header["V47WINDX"])
301  bin = float(hdu[0].header["V47BINX"])
302 
303  hdu.close()
304 
305  rotime = 0
306  if(bin == 1):
307  if(window == 1024):
308  if(pxrt == 2500):
309  rotime = 0.283446688
310  if(pxrt == 250):
311  rotime = 1./0.44
312  if(pxrt == 80):
313  rotime = 1./0.144
314  if(window == 512):
315  if(pxrt == 2500):
316  rotime = 1./6.70
317  if(pxrt == 250):
318  rotime = 1./1.487
319  if(pxrt == 80):
320  rotime = 1./0.535
321  if(window == 256):
322  if(pxrt == 80):
323  rotime = 1./1.772
324  if(window == 64):
325  if(pxrt == 2500):
326  rotime = 1./31.480
327  if(window == 32):
328  if(pxrt == 2500):
329  rotime = 1./42.779
330  if(bin == 2):
331  if(window==1024):
332  if(pxrt==80):
333  rotime = 1./0.551
334  if(bin==16):
335  if(window==1024):
336  if(pxrt==80):
337  rotime=1./10.420
338 
339  if(rotime == 0):
340  print 'Read out time not found for %s' % fname
341  return [-1,-1]
342 
343 
344 
345  return [rotime, exptime]
346 
347 def getLoopStat(datadir, aosysdir, force = 0):
348  """
349  Calculates loop status and avg WFE during each image in a directory.
350  Loop status is closed iff the loop was closed during the entire exposure. If exposure time is less than 1 second, the nearest time is used.
351  A check of WFE is first made for inf and spuriously large values, which are interpolated. Avg WFE is calculated as the average of the 1 second averages recorded in the system logs.
352  If exptime is shorter than 1 second, the nearest WFE entry is used. Std dev of WFE is calculated as the average of the 1 second averages recorded in the system logs.
353 
354  Writes a data file to the datadir with the wfe, and a file with the second by second system status.
355 
356  input:
357  datadir = directory where the images are located
358  aosysdir = directory where the ao system logs are located
359  force = perform update even if the file has already been updated.
360 
361  output:
362  list of [full_path, file_name, loop_status, wfe_avg, std_avg]
363  """
364 
365  #First load a single file just to see if this directory has been processed.
366  datalistone = getOneFileTime(datadir)
367 
368  if len(datalistone) < 1:
369  print 'No VisAO fits files found in %s' % datadir
370  return (['-1', '-1', '-1', '-1', '-1'])
371 
372  if checkLoopStat(datalistone) == 1 and force == 0:
373  return (['-2', '-2', '-2', '-2', '-2'])
374 
375 
376  datalist = getFileTimes(datadir)
377 # if checkLoopStat(datalist) == 1 and force == 0:
378 # return (['-2', '-2', '-2', '-2', '-2'])
379 
380  loglist = getFileTimes(aosysdir,'aosys_', '*.txt')
381 
382  #print datalist
383 
384  if len(datalist) < 1:
385  print 'No visao fits files found'
386  return (['-1', '-1', '-1', '-1', '-1'])
387 
388  print 'Getting system logs:'
389  print ' Directory: %s' % datadir
390  print ' Start time: %f' % datalist[0][2]
391  print ' End time: %f' % datalist[len(datalist)-1][2]
392  print ' Elapsed time: %f' % (float(datalist[len(datalist)-1][2]) - float(datalist[0][2]))
393 
394  #get the log entries that apply to this data-set
395  #todo: make gLE just use times, and calculate ftime1 here from datalist[0][3] and [4]
396  rawentries = getLogEntries(datalist[0][2]-datalist[0][3]-datalist[0][4], datalist[len(datalist)-1][2], loglist)
397 
398  #make a deep copy so we can compare interpolated vs raw values.
399  entries = copy.deepcopy(rawentries)
400 
401  print 'Starting WFE interpolation'
402  #check for and interpolate bad WFE values
403  for i in range(1,len(entries)-1):
404  if entries[i][2] == 'inf' or float(entries[i][2]) > 1e5:
405  n = 0
406  q = 1
407  while n == 0:
408  if(entries[i-q][0] != 'inf' and float(entries[i-q][2] < 1e5)):
409  n = q
410  q = q+1
411  if i-q < 0:
412  n=i
413 
414  t0 = float(entries[i-n][0])
415  wfe0 = float(entries[i-n][2])
416  std0 = float(entries[i-n][3])
417 
418  t1 = float(entries[i][0])
419 
420  if(entries[i+1][0] != 'inf' and float(entries[i+1][2] < 1e5)):
421  n = 1
422  else:
423  n=2
424 
425  n = 0
426  q = 1
427  while n == 0:
428  if(entries[i+q][0] != 'inf' and float(entries[i+q][2] < 1e5)):
429  n = q
430  q = q+1
431  if i+q > len(entries)-1:
432  n = len(entries)-1 - i
433 
434 
435  t2 = float(entries[i+n][0])
436  wfe2 = float(entries[i+n][2])
437  std2 = float(entries[i+n][2])
438 
439  wfe1 = wfe0 + (wfe2-wfe0)/(t2-t0)*(t1-t0)
440  std1 = std0 + (std2-std0)/(t2-t0)*(t1-t0)
441 
442  entries[i][2] = str(wfe1)
443  entries[i][3] = str(std1)
444 
445  #write the wfe time series
446  f = file(datadir+'/wfe.txt', 'w')
447  for i in range(len(rawentries)):
448  f.write("%s %s %s %s %s %s \n" % (rawentries[i][0], rawentries[i][1], rawentries[i][2], rawentries[i][3], entries[i][2], entries[i][3]))
449 
450  #now get the system status for each exposure
451  statent = list()
452  for i in range(len(datalist)):
453 
454  joetime = getJoeTimes(datalist[i][0])
455 
456  etimes = getVisAOExpTimes(datalist[i][2], joetime[0], joetime[1])
457 
458  if len(entries) <= 0:
459  print 'No AOSYS entries'
460  return (['-1', '-1', '-1', '-1', '-1'])
461 
462  #print "%s %f %f %s" %(entries[0][0], etimes[0], etimes[1], entries[len(entries)-1][0])
463  j = 0
464  while(float(entries[j][0]) <= etimes[0]):
465  j=j+1
466  if j >= len(entries): break
467 
468  #j = j
469 
470  loopst = 1
471  wfeavg = 0.
472  stdavg = 0.
473 
474  k = j
475 
476  if k < len(entries):
477 
478  while(float(entries[k][0]) <= etimes[1] and k < len(entries)-1):
479 
480  if(int(entries[k][1]) != 1): loopst = 0
481  wfeavg = wfeavg + float(entries[k][2])
482  stdavg = stdavg + float(entries[k][3])
483  k = k + 1
484  if k >= len(entries): break
485 
486  if(k < len(entries)-1):
487 
488  if(k == j):
489  if(int(entries[j][1]) != 1 or int(entries[j+1][1]) != 1): loopst = 0
490  wfeavg = 0.5*(float(entries[j][2]) + float(entries[j][2]))
491  if(j < len(entries)-1):
492  stdavg = 0.5*(float(entries[j][3]) + float(entries[j+1][3]))
493  else:
494  stdavg = 0.
495  else:
496  wfeavg = wfeavg/(k-j)
497  stdavg = stdavg/(k-j)
498 
499  statent.append([datalist[i][0], datalist[i][1], loopst, wfeavg, stdavg])
500 
501  f = file(datadir+'/status.txt', 'w')
502  for i in range(len(statent)):
503  f.write("%s %i %f %f\n" % (statent[i][1], statent[i][2], statent[i][3], statent[i][4]))
504 
505  return statent
506 
507 def updateLoopStat(statlist):
508 
509  #print statlist
510 
511  for im in statlist:
512  #print im[0]
513  hdu = pyfits.open(im[0], 'update')
514 
515  if im[2] == 1:
516  hdu[0].header.update('AOLOOPST', "CLOSED",comment="AO loop status during exposure")
517  #hdu[0].header.comments['AOLOOPST'] =
518  else:
519  hdu[0].header.update('AOLOOPST', "OPEN",comment="AO loop status during exposure")
520 
521  hdu[0].header.update('AVGWFE', im[3],comment='Avg WFE (nm rms phase)')
522  hdu[0].header.update('STDWFE', im[4],comment='Std Dev of WFE (nm rms phase)')
523 
524  hdu[0].header.update('HISTORY', 'VisAO system status update applied %s' % time.asctime(time.localtime(time.time())))
525 
526  hdu.close()
527 
528 def checkLoopStat(statlist):
529  #Check if this list has already been processed
530 
531  if len(statlist) <= 0:
532  return 0
533 
534  im = statlist[0]
535  hdu = pyfits.open(im[0], 'update')
536 
537  if hdu[0].header.get('AOLOOPST') != 'NOT PROCESSED':
538  return 1
539 
540  return 0
541 
542 
543 def getGoodIms(dir, movedir):
544  """
545  run this after headers updated
546  """
547  if(dir[len(dir)-1] != '/'):
548  dir = dir + '/'
549 
550  if(movedir[len(movedir)-1] != '/'):
551  movedir = movedir + '/'
552 
553  fnames = glob.glob(dir+"V47_*.fits")
554 
555  good = list()
556 
557  nopen = 0
558  nsci = 0
559  ndrk = 0
560 
561  for dfn in fnames:
562  hdu = pyfits.open(dfn)
563  fn = dfn[len(dir):len(dfn)]
564 
565  if (hdu[0].header["AOLOOPST"] == 'OPEN' and hdu[0].header["VIMTYPE"] == 'SCIENCE'):
566  nopen = nopen + 1
567 
568 
569  if (hdu[0].header["AOLOOPST"] == 'CLOSED' and hdu[0].header["VIMTYPE"] == 'SCIENCE'):
570  nsci = nsci + 1
571  good.append([dfn, fn])
572 
573  if (hdu[0].header["VIMTYPE"] == 'DARK'):
574  ndrk = ndrk + 1
575  good.append([dfn, fn])
576 
577  print nopen
578  print nsci
579  print ndrk
580 
581 
582  for f in good:
583  newf = movedir + f[1]
584  shutil.copyfile(f[0], newf)
585 
586 
587 def getDarks(dir, movedir):
588  """
589  run this after headers updated
590  """
591  if(dir[len(dir)-1] != '/'):
592  dir = dir + '/'
593 
594  if(movedir[len(movedir)-1] != '/'):
595  movedir = movedir + '/'
596 
597  fnames = glob.glob(dir+"V47_*.fits")
598 
599  good = list()
600 
601  nopen = 0
602  nsci = 0
603  ndrk = 0
604 
605  for dfn in fnames:
606  hdu = pyfits.open(dfn)
607  fn = dfn[len(dir):len(dfn)]
608 
609  if (hdu[0].header["VIMTYPE"] == 'DARK'):
610  ndrk = ndrk + 1
611  good.append([dfn, fn])
612 
613  print nopen
614  print nsci
615  print ndrk
616 
617 
618  for f in good:
619  newf = movedir + f[1]
620  shutil.copyfile(f[0], newf)
621 
622 def sdssFilters():
623  """
624  Returns a list of the standard VisAO filters
625  """
626  filters=["r'", "i'", "z'", "Ys"]
627  return filters
628 
629 def xDither(scale = 1.0):
630  """
631  Returns a list of x offsets for a gimbal dither. Units are mm.
632 
633  scale = total width of the pattern in mm
634  """
635  dx = [-.5*scale, 1.*scale, 0.*scale, -1.*scale]
636 
637  return dx
638 
639 def yDither(scale = 1.0):
640  """
641  Returns a list of y offsets for a gimbal dither. Units are mm.
642 
643  scale = total width of the pattern in mm
644  """
645  dy = [.5*scale, 0.*scale, -1.*scale, 0.*scale]
646 
647  return dy
648 
649 def visao_say(txt):
650  cmd = "echo vicki %s | nc zorro.lco.cl 50000" % (txt)
651  os.system(cmd)
652 
653 def visao_script_complete():
654  visao_say("viz ay oh script complete")
655 
656 def visao_script_error():
657  visao_say("viz ay oh script airor")
658 
659 def visao_script_stopped():
660  visao_say("viz ay oh script stopped")
661 
662 
663