spectractor.midworker

1 2 '''Spectractor submodule containing classes which interact with ESO MIDAS. 3 Interaction with ESO MIDAS was made through pymidas module. 4 With pyMidas we calculate heliocentric correction and 5 with pyephem we calculate precession. 6 Midworker class is the main part submodule. 7 8 @author: Dmitry Nasonov 9 ''' 10 11 import os, sys, shutil 12 from datetime import datetime 13 from math import pi 14 15 import numpy as np 16 17 try: 18 from pymidas import midas 19 except ImportError, err: 20 print >> sys.stderr, str(err) + ": We use ESO MIDAS and PyMidas for, eg. ", 21 print "heliocentric RV correction. Get it here: http://www.eso.org/esomidas" 22 sys.exit(True) 23 24 try: 25 from coords.astrodate import utc2jd, jd2jyear 26 except ImportError, err: 27 print >> sys.stderr, str(err) + ": coords module required, ", 28 print >> sys.stderr, "download and install it from ", 29 print >> sys.stderr, "http://www.scipy.org/AstroLib/ ", 30 print >> sys.stderr, "or https://www.stsci.edu/trac/ssb/astrolib" 31 # Another link: http://stsdas.stsci.edu/astrolib/coords_api/ 32 sys.exit(True) 33 34 try: 35 from ephem import Equatorial 36 #import ephem._libastro as libastro 37 except ImportError, err: 38 print >> sys.stderr, str(err) + ": pyephem module required ", 39 print >> sys.stderr, "for coordinates precessing." 40 41 try: 42 import pyfits 43 except ImportError, err: 44 print >> sys.stderr, str(err), "PyFITS is required for read and write FITS" 45 46 from serve import Logger 47 48 spacer = lambda *args: " %s " % ' '.join(map(str, args)) 49 commer = lambda arg: ",".join(map(str, arg)) 50 convcoord = lambda cr, rest: (int(cr), int(divmod(rest*60, 1)[0]), 51 round(divmod(rest*60, 1)[1]*60, 3)) 52 53

54 -class HelVelCor:

55 """Calculate heliocentric RV correction for number of observation times. 56 @param params: list of observation parameters: date, time, coordinates. 57 See docstrings of methods of this class for details 58 @param dtype: data type of resulting array, must be dict or list 59 @param obscoord: Observatory coordinates. Default is for 6-m telescope 60 (BTA), SAO RAS 61 """

62 - def __init__(self, params, dtype, obscoord=None, verbose=False):

63 self.params = params 64 self.verbose = verbose 65 if not obscoord: 66 obscoord = ((41, 26, 30), (43, 39, 12)) # 41.4414 43.6467 67 self.obscoord = spacer(commer(obscoord[0]), commer(obscoord[1])) 68 self.dtype = dtype 69 self.fdata = dtype() 70 self.intcoord = lambda h, m, s: (int(h), int(m), float(s)) 71 self.cmpbary = lambda date, time: ' '.join(map(','.join, ( 72 map(str, date), map(str, time))))

73

74 - def precess_midas(self, date, time, objcoords):

75 """Wrapper for comp/prec MIDAS routine. 76 @param date: tuple of year, month and day 77 @param time: tuple of hours and minutes 78 @return: ra and dec as tuples 79 """ 80 minuhour = (time[0]/24.) + (time[1]/1440.) 81 compspec_str = ",".join(map(str, date)) + '.' + \ 82 str(int(minuhour*10000)) 83 ra, dec, epoch = objcoords 84 coords_str = spacer(commer(ra), commer(dec), str(epoch)) 85 cmpr = midas.computePrec(coords_str, compspec_str) 86 ra = self.intcoord(cmpr[6], cmpr[8], cmpr[10]) 87 dec = self.intcoord(cmpr[13], cmpr[15], cmpr[17]) 88 #jyear_epoch = float(cmpr[3]) 89 return ra, dec

90

91 - def compbary(self, date, time, ra, dec):

92 """Wrapper for comp/bary MIDAS routine. 93 @param date: tuple of year, month and day 94 @param time: tuple of hours and minutes 95 @param ra: tuple of coords 96 @param dec: tuple of coords 97 @return: heliocentric and barycentric radial velocity corrections 98 """ 99 compbary_str = self.cmpbary(date, time) 100 coords_str = spacer(','.join(map(str, ra)), ','.join(map(str, dec))) 101 cmbr = midas.computeBary(compbary_str, coords_str, self.obscoord) 102 hcor = round(float(cmbr[20]), 4) 103 bcor = round(float(cmbr[14]), 4) 104 return hcor, bcor

105

106 - def filler(self, ID, jd, jyr, ra, dec, brcor, hlcor):

107 """Dummy function for filling resulting array""" 108 if self.dtype is dict: 109 self.fdata[ID] = (jyr, round(jd-2400000.5, 5), (ra, dec), 110 brcor, hlcor) 111 elif self.dtype is list: 112 self.fdata.append((ID, jyr, round(jd-2400000.5, 5), (ra, dec), 113 brcor, hlcor))

114

115 - def getdatime(self, spec):

116 """Convert datetime of observetion to JD. Pay attention to indices! 117 """ 118 if type(spec[0]) is not type(datetime(1970, 1, 1)): 119 date, time = spec[:3], spec[3:5] 120 jd = utc2jd(datetime(*spec[:5])) 121 else: 122 date, time = spec[0].timetuple()[:5] 123 jd = utc2jd(spec[0]) 124 return date, time, jd

125

126 - def veltimes(self, objcoords, midprecess=False):

127 """Compute heliocentric correction for the set of observation times. 128 @param objcoords: Object coordinates with its epoch as tuple: 129 (ra, dec, epoch). ra and dec are also tuples. 130 @param midprecess: Precess coordinates using ESO MIDAS comp/bary 131 routine, use pyephem instead 132 @return: list of (spectra) parameters, including heliocentric and 133 barycentric radial velocity correction 134 """ 135 for ID, spec in sorted(self.params.items()): 136 date, time, jd = self.getdatime(spec) 137 jyr = round(jd2jyear(jd), 5) 138 if midprecess or abs(objcoords[2]-2000) > .01: 139 # Precess coords to jyear_epoch using MIDAS COMP/PREC 140 ra, dec = self.precess_midas(date, time, objcoords) 141 elif not midprecess: 142 # Precess coords using pyephem 143 ra, dec = precess(jd-2400000.5, objcoords, verbose=self.verbose) 144 if self.verbose: 145 print "Precessed coordinates:", ra, dec 146 helrvcor, barrvcor = self.compbary(date, time, ra, dec) 147 self.filler(ID, jd, jyr, ra, dec, barrvcor, helrvcor) 148 return self.fdata

149

150 - def velcoords(self, midprecess=False):

151 """Compute heliocentric correction for the set of coordinates. 152 """ 153 for ID, spec in sorted(self.params.items()): 154 date, time, jd = self.getdatime(spec) 155 jyr = round(jd2jyear(jd), 5) 156 objcoords = spec[5:8] 157 if not objcoords[2]: 158 objcoords[2] = jyr 159 if midprecess or abs(objcoords[2]-2000) > .01: 160 # Precess coords to jyear_epoch using MIDAS COMP/PREC 161 ra, dec = self.precess_midas(date, time, objcoords) 162 elif not midprecess: 163 # Precess coords using pyephem 164 ra, dec = precess(jd-2400000.5, objcoords, verbose=self.verbose) 165 if self.verbose: 166 print "Precessed coordinates:", ra, dec 167 helrvcor, barrvcor = self.compbary(date, time, ra, dec) 168 self.filler(ID, jd, jyr, ra, dec, barrvcor, helrvcor) 169 return self.fdata

170 171

172 -def precess(mjd, objcoords, verbose=False):

173 """Precess coords using pyephem from J2000 to mjd. 174 @param mjd: Modified julian date: JD - 2400000.5 175 @param objcoords: Object coordinates with its epoch as tuple 176 @return: ra, dec as tuples 177 """ 178 # Dublin Julian Day = JD - 2415020.0 179 #J2000 = 2451545.0 - 2415020.0 <-- it is from libastro sources 180 # Maybe we should try direct call to libastro: 181 #libastro.precess(2000, 2007.794, 1.386, 0.282) 182 #2415020.0 - 2400000.5 = 15019.5 183 mjd = mjd - 15019.5 184 if verbose: 185 print "We get the coords", objcoords 186 # coords as strings, separated by ":" 187 dv = lambda arg: ":".join(map(str, arg)) 188 coord = Equatorial(dv(objcoords[0]), dv(objcoords[1])) 189 precessed_coords = Equatorial(coord, epoch=mjd) 190 ra, dec = precessed_coords.to_radec() 191 # Get tuple for each coordinate 192 ra, dec = convcoord(*divmod(ra*12/pi, 1)), convcoord(*divmod(dec*180/pi, 1)) 193 if verbose: 194 print "Precess to", mjd, ra, dec, "with J2000 coords", objcoords 195 return ra, dec

196 197

198 -class MidWorker:

199 """Class with simple wrappers for some ESO MIDAS routines 200 """

201 - def __init__(self, usedisp=False, verbose=False):

202 self.verbose = bool(verbose) 203 self.usedisp = usedisp 204 self.ext = 'bdf' 205 self.med = 35 206 self.avermed = lambda m: ",".join(('median', str(m), str(m), 'data'))

207

208 - def plot(self, bnam, plt="row", cut=2010, over=False):

209 """Wrapper for Midas functions plot/row, plot/colu""" 210 if over: 211 over = "over" 212 else: 213 over = "" 214 if self.verbose: 215 print "Plot", over, plt, bnam, "cut on pixel", cut 216 cmd = " ".join((over+"plot/"+plt, bnam, "@"+str(cut/2))) 217 midas.do(cmd)

218

219 - def readesc(self, fname, desc="npix"):

220 """Wrapper for read/desc Midas procedure""" 221 nx, ny = midas.readDesc(fname, desc)[-2:] 222 return int(nx), int(ny)

223

224 - def loima(self, img, chan=0):

225 """Load image with autocuts""" 226 if not (os.path.isfile(img) or os.path.isfile(img+'.'+self.ext)): 227 return None, None 228 nx, ny = self.readesc(img) 229 dx, dy = map(int, midas.writeOut('{ididev(2)} {ididev(3)}')) 230 if nx > dx: 231 i = round(- (nx/(dx - 10.)) - 1, 2) 232 if ny > dy: 233 j = round(- (ny/(dy - 10.)) - 1, 2) 234 scale = ",".join((str(i), str(j), 'a')) 235 if self.verbose: 236 print "Scale to", scale 237 midas.loadImag(img, str(chan), 'scale='+scale, 'center=c,c') 238 midas.writeKeyw("log/i/4/1", "2") 239 #midas.putKeyword("log", 2) 240 midas.statistImag(img) 241 midas.do('@a autocuts ' +img+ ' SIGMA') 242 midas.writeKeyw("log/i/4/1", "0") 243 return i, j

244

245 - def savebdf(self, name, bdf_pth, wridesc=True):

246 """Save a MIDAS bdf image, using existing FITS file 247 Also write descriptors, load it on display 248 """ 249 midas.indiskFits(name, bdf_pth) 250 if wridesc: 251 midas.writeDesc(bdf_pth, 'start/d/1/2 1.,1.') 252 midas.writeDesc(bdf_pth, 'step/d/1/2 1.,1.') 253 if self.usedisp: 254 self.loima(bdf_pth)

255

256 - def crebias(self, files, biasname="bias", med=None, delfits=False):

257 """Create bias frame by averaging a set of images. 258 @param med: value of low and high interval in median averaging 259 """ 260 if not med: 261 med = self.med 262 catname = biasname+".cat" 263 midas.createIcat(catname, "NULL") 264 for file_pth in files: 265 midas.addIcat(catname, file_pth) 266 midas.averageImag(biasname, "=", catname, "? +", self.avermed(med)) 267 self.pltgra(biasname, "row", 1000, over=False) 268 midas.outdiskFits(biasname, biasname+'.fits')

269

270 - def creflat(self, files, flatname="flat", delcat=False, averopt=""):

271 """Create flat field average image 272 @param delcat: delete created cat, if icat is present 273 """ 274 catname = flatname+".cat" 275 midas.createIcat(catname, "NULL") 276 for file_pth in files: 277 midas.addIcat(catname, file_pth) 278 midas.averageImag(flatname, "=", catname, "?", averopt) 279 self.pltgra(flatname, "col", 1000, over=False)

280

281 - def filtcosm(self, specname, ns=2, cosm="cosm"):

282 """Wrapper for filter/cosm Midas routine, which removes cosmic 283 ray events from CCD image and replace them by a local median value. 284 @param specname: input image name. Resulting image will be with the 285 same name and the 'f' postfix 286 @param ns: threshold for the detection of cosmic rays 287 param cosm: frame containing mask for the giving image 288 """ 289 #if specname.endswith(".bdf"): 290 resnam = specname.split(".")[0]+"f" 291 params = ",".join(map(str, (0, 2, 8, ns, 2))) 292 if self.usedisp: 293 i, j = self.loima(specname) 294 midas.filterCosm(specname, resnam, params, cosm) 295 if self.usedisp: 296 scale = ",".join((str(i), str(j), 'max')) 297 midas.loadImag(cosm, "scale="+scale, "cuts=0,1")

298

299 - def cycle(self, beg, end, plt, fcat=False, pfix=""):

300 k = True 301 if fcat: 302 midas.createIcat(fcat, "NULL") 303 else: 304 namlst = [] 305 for bnum in xrange(beg, end+1): 306 bnam = 'l'+str(bnum).zfill(3) + pfix 307 if os.path.isfile(bnam+'.'+self.ext): 308 if k: 309 if self.verbose: 310 print "Add first", bnam, "image" 311 nx, ny = self.readesc(bnam) 312 cut = ny if plt == "row" else nx 313 self.plot(bnam, plt, cut) 314 firstnam, k = bnam, False 315 else: 316 if self.verbose: 317 print "Add", bnam, "image" 318 self.plot(bnam, plt, cut, over=True) 319 if fcat: 320 midas.addIcat(fcat, bnam) 321 else: 322 namlst.append(bnam) 323 if not fcat: 324 namlst = ",".join(namlst) 325 else: 326 namlst = fcat 327 return firstnam, namlst, cut

328

329 - def pltgra(self, resnam, plt, cut, over=True, midlo=False):

330 midas.setGrap("color=2") 331 self.plot(resnam, plt, cut, over=over) 332 midas.setGrap("color=1") 333 if midlo: 334 midas.do("@@ lo_ima " + resnam) 335 elif self.usedisp: 336 self.loima(resnam)

337

338 - def crebiasref(self, beg, end, resnam='bias', icat=False, plt="row", 339 med=None):

340 """Create bias frame from set of biases. 341 @param resnam: resulting frame name 342 @param med: value of low and high interval in median averaging 343 @param plt: plotting mode, could be row or col. 344 """ 345 firstnam, namlst, cut = self.cycle(beg, end, plt, fcat=icat) 346 if not med: 347 med = self.med 348 midas.averageImag(resnam, "=", namlst, "? +", self.avermed(med)) 349 self.pltgra(resnam, plt, cut)

350

351 - def creflatref(self, beg, end, resnam=None, icat='flat.cat', plt="col", 352 pfix="d", averopt="", delcat=False):

353 """Create flat field average image 354 @param delcat: delete created cat, if icat is present 355 @param pfix: postfix of image filename, the "d" value means 356 dark substracted 357 """ 358 firstnam, namlst, cut = self.cycle(beg, end, plt, fcat=icat, pfix=pfix) 359 if not resnam: 360 resnam = firstnam.replace('d', 'ad') 361 midas.averageImag(resnam, "=", namlst, "?", averopt) 362 self.pltgra(resnam, plt, cut) 363 if delcat and icat: 364 os.remove(icat)

365

366 - def subdark(self, beg, end, biasname="bias", pfix="d"):

367 """Substract bias from images""" 368 for bnum in xrange(beg, end+1): 369 bnam = 'l'+str(bnum).zfill(3) 370 if os.path.isfile(bnam+'.'+self.ext): 371 nameo = bnam + pfix 372 midas.computeImag(nameo, "=", bnam, "-", biasname) 373 midas.statistImag(nameo) 374 if self.usedisp: 375 self.loima(nameo) 376 midas.copyDd(bnam, "*,3", nameo) 377 os.remove(bnam+'.'+self.ext)

378 379

380 -class Preparer(MidWorker, Logger):

381 """Class for image processing with MIDAS, against MIDAS sripting language =) 382 We use MIDAS routines through pyMidas, but replace it with numpy and pyfits, 383 if it is possible. 384 385 This code also heavily depends of some specific characteristics of SAO RAS 386 spectral archive, mostly on spectra taken on Nesmyth Echelle Spectrograph of 387 6-meter telescope. 388 It is not pipeline package, but a set of non-interactive methods for packet 389 processing. It contains methods for averaging images, substracting bias 390 frame etc. 391 @param dest_dir: destination dir for moving processed fits files. Could be 392 omitted. 393 @param usedisp: use Midas display for loading images using Midas routines 394 """

395 - def __init__(self, dest_dir="arch", usedisp=False, verbose=False):

396 self.dest_dir = dest_dir 397 self.verbose = bool(verbose) 398 self.usedisp = usedisp 399 self.fitsfixmode = {True: 'fix', False: 'silentfix'}[self.verbose] 400 self.sep, self.postfix, self.ext = '_', '.', 'fits' 401 #self.imgcut_dct = {2052: (4, 1992, 4, 2039), # NES 2052x2052 chip 402 self.imgcut_dct = {2052: (3, 2048, 4, 2039), # NES 2052x2052 chip 403 # Resulting shape after cut for NES images with cuts (4, 1992, 4, 2039) 404 # will be (2035, 1988). NES with full use: (3, 2048, 0, 2046), 405 # shape will be (2045, 2046) 406 2068: (0, 2048, 2, 2048), # LPR 2068x2072 chip 407 1050: (5, 1040, 2, 1157) # PFES/Lynx 1050x1170 chip 408 } 409 self.nesbadrows = [(338, 337, 339), (568, 567, 569), 410 (336, 336, 335), (335, 336), (1044, 1044, 1043, 1042), 411 (1043, 1044), (1042, 1044), (433, 433, 434, 435), 412 (434, 433), (435, 433), (185, 185, 184), (184, 185), (1903, 1903, 1902), 413 (1902, 1903), (1229, 1229, 1228), (1228, 1229), (47, 47, 48), (48, 47)] 414 self.rowdoc = lambda d, xs: sum([d[:, 1992-r] for r in xs]) / len(xs) 415 self.mw = MidWorker(usedisp=usedisp, verbose=verbose) 416 self.Lg = Logger(fitsfixmode=self.fitsfixmode, verbose=verbose)

417

418 - def prepare(self, files, log=True, crop=True, rot=True, badrow=False, 419 substbias=False, cleanhead=True, filmove=True, process=True):

420 """Method for packet preparing and/or logging of raw images, 421 mostly taken on spectrographs of 6-meter telescope. 422 @param files: list of raw image parameters with its paths 423 @param log: create log, if True. if "logonly", create log 424 without image processing 425 @param crop: crop each image depending on its shape, which is 426 determine used CCD chip 427 @param rot: rotate image clockwise 428 @param badrow: mask bad rows on NES CCD chip 429 @param substbias: substract bias frame, if the filename 430 bias image is given 431 @param cleanhead: remove empty strings in fits headers 432 @param filmove: move or not source files (raw images) to 433 self.arch_dir, after they are have been processed 434 @param process: process images, if True. Otherwise all 435 of the processing procedures (crop, rot, badrow, substbias) 436 will be omitted 437 """ 438 if substbias and process and type(substbias) is str: 439 self.bias = pyfits.getdata(substbias) 440 self.substbias = True 441 else: 442 self.substbias = False 443 self.crop = crop 444 self.rot = rot 445 self.badrow = badrow 446 if files and filmove and not os.path.isdir(self.dest_dir): 447 os.mkdir(self.dest_dir) 448 for (file_pth, num, date, fflag, keymode) in sorted(files): 449 if self.verbose: 450 print pyfits.info(file_pth) 451 print "Open file", file_pth 452 curfts = pyfits.open(file_pth) #mode="update" 453 if cleanhead: 454 # Clean header - delete empty strings: 455 del curfts[0].header[''] 456 if log: 457 self.Lg.logger(curfts, num, date, file_pth, keymode=keymode) 458 if log == "logonly": 459 continue 460 if file_pth in self.Lg.biases: 461 if filmove: 462 shutil.move(file_pth, self.dest_dir) 463 continue 464 elif file_pth in self.Lg.flats: 465 newname = self.Lg.flats[file_pth] 466 hilim = 9000 467 elif file_pth in self.Lg.calib: 468 newname = self.Lg.calib[file_pth] 469 hilim = 5000 470 else: 471 date = "".join((str(date[0]), str(date[1]).zfill(2), 472 str(date[2]).zfill(2))) 473 newname = os.path.join(self.dest_dir, "".join((date, self.sep, 474 num, self.postfix, self.ext))) 475 hilim = 2600 476 if process: 477 self.processdata(curfts, hilim=hilim) 478 self.savefits(curfts, newname) 479 if filmove: 480 shutil.move(file_pth, self.dest_dir) 481 if file_pth in self.Lg.flats: 482 continue 483 elif file_pth in self.Lg.calib: 484 self.mw.filtcosm(newname, cosm="c"+num) 485 continue 486 self.mw.savebdf(newname, "l"+num+"d.bdf")

487

488 - def crebias(self, biasname="bias", files=None, filmove=True, log=True):

489 """Average bias frames with crebias Midworker method""" 490 if files: 491 self.prepare(files, process=False, cleanhead=False, 492 log=log, filmove=filmove) 493 elif not self.Lg.biases: 494 print "No biases, return" 495 return 496 for file_pth, newname in self.Lg.biases.items(): 497 if not os.path.isfile(newname): 498 shutil.copy2(file_pth, newname) 499 self.mw.crebias(self.Lg.biases.values(), biasname=biasname)

500

501 - def creflat(self, flatname="flat"):

502 """Average flat field frames with creflat Midworker method. 503 After averaging make a filter/cosm procedure 504 """ 505 if self.Lg.flats: 506 self.mw.creflat(self.Lg.flats.values(), flatname=flatname) 507 self.mw.filtcosm(flatname)

508

509 - def processdata(self, curfts, hilim=600):

510 """Process image data: substract bias, crop, rotate, mask bad rows""" 511 data = curfts[0].data 512 dshape = data.shape 513 if self.substbias: 514 data -= self.bias 515 if self.badrow: 516 self.nesbadrow(data, hilim=hilim) 517 if self.crop and dshape[0] in self.imgcut_dct: 518 x1, x2, y1, y2 = self.imgcut_dct[dshape[0]] 519 #midas.extractImag({name} = {namef}[@{x1},@{y1}:@{x2},@{y2}]) 520 data = data[y1:y2, x1:x2] 521 if self.rot and dshape[0] >= 2052: 522 data = np.rot90(data) 523 #if self.flip: 524 ## flip for Lynx/Pfes 525 #data = np.flipud(data) 526 # Write result 527 curfts[0].data = data 528 curfts[0].scale('float32', 'old') # int16

529

530 - def savefits(self, curfts, newname):

531 """Save current fits""" 532 if os.path.isfile(newname): 533 bcknam = newname.replace("."+self.ext, "-backup."+self.ext) 534 shutil.move(newname, bcknam) 535 try: 536 curfts.writeto(newname) 537 except pyfits.core.VerifyError: 538 curfts.verify(self.fitsfixmode) 539 curfts.writeto(newname) 540 curfts.close()

541

542 - def nesbadrow(self, d, lowcut=-25, hilim=600):

543 """Mask bad pixels in raw images taken with Uppsala NES chip. 544 """ 545 for row in self.nesbadrows: 546 d[:, 2052-row[0]-60] = self.rowdoc(d, row[1:]) 547 d[:, 1992] = (d[:, 1991] + d[:, 1993]) / 2 548 d[:, 1997] = (d[:, 1996] + d[:, 1998]) / 2 549 d[:, 2005] = (d[:, 2004] + d[:, 2006]) / 2 550 d[:, 2008] = (d[:, 2007] + d[:, 2009]) / 2 551 d[:, 2010] = (d[:, 2009] + d[:, 2011]) / 2 552 d[:, 2014] = (d[:, 2013] + d[:, 2015]) / 2 553 d[:, 2027] = (d[:, 2026] + d[:, 2028]) / 2 554 for i in xrange(2036, 2038): 555 d[i, :22] = (d[i-1, :22] + d[i-2, :22])/2 556 d[i, 1337:1569] = (d[i-1, 1337:1569] + d[i-2, 1337:1569])/2 557 for i in xrange(2038, 2040): 558 d[i, :24] = (d[i-1, :24] + d[i-2, :24])/2 559 d[i, 1337:1569] = (d[i-1, 1337:1569] + d[i-2, 1337:1569])/2 560 for i in xrange(2040, 2042): 561 d[i, :26] = d[i-1, :26] 562 d[i, 1125:1692] = (d[i-1, 1125:1692] + d[i-2, 1125:1692])/2 563 for i in xrange(2042, 2044): 564 d[i, :34] = d[i-1, :34] 565 #d[i, 1087:27] = d[i-1, 1087:27] 566 #for i in xrange(2044, 2046): 567 #d[i, :2052-2005] = d[i-1, :2052-2005] 568 #d[i, 2052-1130:2052-15] = d[i-1, 2052-1130:2052-15] 569 #zonlim = (0, 250, 1979) #[31:96, 1988:2040] 570 zonlim = (1, 118, 1988) 571 zone = d[zonlim[0]:zonlim[1], zonlim[2]:] 572 med = min(np.median(zone), hilim) 573 hilim = max(np.mean(d), hilim) 574 d[zonlim[0]:zonlim[1], zonlim[2]:][np.where(zone > hilim)] = med 575 #d[:280, 2004:2006][np.where(d[:280, 2004:2006] > hilim)] = med 576 d[:280, 2004:2011][np.where(d[:280, 2004:2011] > hilim)] = med 577 d[np.where(d < lowcut)] = lowcut

578

579 - def savelog(self, filenam="night", wdir=""):

580 """Save observation log to file.""" 581 if os.path.isfile(os.path.join(wdir, filenam + ".log")): 582 shutil.move(os.path.join(wdir, filenam + ".log"), 583 os.path.join(wdir, filenam + ".old")) 584 logfile = open(os.path.join(wdir, filenam+".log"), "w") 585 log = self.Lg.logtostr() 586 for line in log: 587 print >> logfile, line 588 print >> logfile 589 log = self.Lg.journaltostr() 590 for line in log: 591 print >> logfile, line 592 logfile.close()

593

Source Code for Module spectractor.midworker