nictools.saaclean

1 ## Automatically adapted for numpy.numarray Feb 05, 2007 by 2 3 """ 4 saaclean: Module for estimating and removing persistent CR signal due to a prior 5 SAA passage. 6 7 Usage: Normally used via the STSDAS task saaclean in the nicmos package. 8 To use as pure python, create a params object to override any of 9 the default parameters if desired, then invoke clean: 10 >>> mypars=saaclean.params(thresh=0.23) 11 >>> saaclean.clean('inputfile.fits','outputfile.fits',pars=mypars) 12 13 For more information: 14 Additional user information, including parameter definitions and more 15 examples, can be found in the help file for the STSDAS saaclean task, 16 located in nicmos$doc/saaclean.hlp. 17 18 The algorithm and IDL prototype are described in the NICMOS 19 ISR 2003-009, by Bergeron and Dickinson, available through the NICMOS 20 webpage. 21 22 Dependencies: 23 numpy 1.0.2.dev3534 or higher 24 pyfits v1.1b4 or higher 25 imagestats v1.3 or higher 26 27 """ 28 29 from __future__ import division 30 31 __version__="1.3" 32 __vdate__="2009-12-14" 33 34 # The above text is duplicated in the __init__ file for the package, since 35 #that's where it shows up for the user. 36 from pytools import numerixenv #Temporary NUMERIX environment check 37 import os 38 import exceptions 39 import numpy as N, pyfits 40 import imagestats #to check the version 41 from imagestats import ImageStats as imstat #pyssg lib 42 from imagestats.histogram1d import histogram1d 43 import SP_LeastSquares as LeastSquares #Excerpt from Hinsen's Scientific Python 44 from numpy.linalg import LinAlgError 45 46 #History: 47 # New version of imagestats: 14 Dec 09, Laidler 48 # - imagestats v1.3 uses "midpt" to duplicate iraf imstats functionality 49 # Enhancements, 20 Jan 06, Laidler 50 # - replaced infile by calcfile and targfile 51 # - allow applying correction to a file other than that on which it 52 # was computed 53 # - added "clobber" parameter to control clobber behavior on all 54 # output files. 55 # Enhancement, 21 Dec 05, Laidler 56 # use header value of ADCGAIN in place of GAINPLOT 57 # add new CRTHRESH, NOISETHRESH, BINSIGFRAC parameters 58 # Enhancement, 21 Dec 05, Laidler 59 # - add chi2 output to header keyword set. This involved changing 60 # the signature of the parabola_min function to return the chi2. 61 # Enhancement, 19 Dec 05, Laidler 62 # - use DQ extension for bad pixels 63 # Bugfix,21 Jan 05, Laidler 64 # - make all paths & filenames more robust via osfn 65 # Bugfix,20 Jan 05, Laidler: 66 # - correct handling of middle column/row in Exposure.pedskyish() 67 # Bugfixes, 3 Aug 04, Laidler: 68 # - make filename construction more robust via os.path.abspath on the directory 69 # - Ensure directory specified in pars.darkpath exists 70 # - Fall back to extension 0 if extension 1 doesn't exist (Exposure class) 71 # - Support noise reduction in high domain only 72 # Initial python implementation: Dec 2003, Laidler 73 # Based on IDL implementation by Bergeron 74 75 #Notes for future improvement: 76 # - possibly make saaper its own class 77 # - the crthreshholding code is kind of tacky 78 # - possibly make filename its own class so it can have a method for nref 79 # instead of using the osfn helper function 80 #........................................................................ 81 #Class definitions 82 #......................................................................... 83

84 -class params:

85 - def __init__(self,scale=0.54,wf1=0.7,wf2=0.3, 86 stepsize=0.008,thresh=None,hirange=0.4,lorange=0.25,dofit=1, 87 crthresh=0.3, noisethresh=1.0, binsigfrac=0.3, 88 readsaaper='False',writesaaper='True',saaperfile='saaper.fits', 89 fitthresh='True',histbinwidth=0.01, nclip=10, 90 clobber='False', 91 flatsaaper='True',flatsaaperfile=None, 92 maskfile=None,darkpath=None,diagfile=None):

93 self.scale=scale 94 self.wf1=wf1 95 self.wf2=wf2 96 self.writesaaper=writesaaper 97 self.readsaaper=readsaaper 98 self.saaperfile=osfn(saaperfile) 99 self.clobber=clobber 100 self.flatsaaper=flatsaaper 101 self.flatsaaperfile=osfn(flatsaaperfile) 102 self.maskfile=osfn(maskfile) 103 self.stepsize=stepsize 104 self.thresh=thresh 105 self.hirange=hirange 106 self.lorange=lorange 107 self.dofit=dofit 108 109 self.fitthresh=fitthresh 110 self.histbinwidth=histbinwidth 111 self.nclip=nclip 112 113 self.crthresh=crthresh 114 self.noisethresh=noisethresh 115 self.binsigfrac=binsigfrac 116 117 self.darkpath=osfn(darkpath) 118 self.diagfile=osfn(diagfile) 119 120 self.appstring=None # Might be needed later.

121

122 -class Domain:

123 """ Stores a list of pixels for a (typically high or low) signal domain""" 124

125 - def __init__(self,name,pixellist,range):

126 self.name=name 127 self.pixlist=pixellist 128 self.range=range 129 self.npix=len(self.pixlist[0])

130 #Because the pixlist is created with a "where" statement, it's a 131 #2 element array, 1 for x & 1 for y. Thus to get the number of 132 #pixels, we need the length in one of the elements. 133 134 135

136 - def striplowerthan(self,factor):

137 """self.pp is defined in Exposure.getscales 138 It contains the (bin, stddev, mode) for the statistical analysis. 139 striplowerthan(factor) examines the stddev column only, and replaces 140 all values of the stddev that are less than factor*the zeroth bin, 141 with the maximum stddev. 142 """ 143 p1=self.pp[1,:] 144 uu=N.where(p1 < factor*p1[0]) 145 if uu[0].size != 0: 146 p1[uu]=p1.max() 147 self.pp[1,:]=p1

148

149 - def getmin(self):

150 ubest=N.where(self.pp[1,:] == self.pp[1,:].min())[0][0] 151 umode=N.where(self.pp[2,:] == self.pp[2,:].min())[0][0] 152 return ubest, umode

153

154 - def writeto(self,filename,clobber=False):

155 if not clobber: 156 if os.path.exists(filename): 157 raise IOError, "%s already exists: aborting\n"%filename 158 #if clobber=True or file does not exist, proceed anyhow 159 f=open(filename,'w') 160 f.write('# '+self.name+'\n') 161 f.write('# Pixels in this domain: '+`len(self.pixlist[0])`+'\n') 162 f.write('# 1 scale factor \n') 163 f.write('# 2 sigma \n') 164 f.write('# 3 mode \n') 165 for i in range(len(self.pp[0])): 166 f.write('%f %f %f\n' % (self.pp[0,i],self.pp[1,i],self.pp[2,i])) 167 f.close()

168 169

170 -class Exposure:

171 """ Stores a collection of keywords and the image data for an exposure. """ 172

173 - def __init__(self,imgfile,nickname=None):

174 175 self.filename=osfn(imgfile) 176 if nickname is None: 177 self.nickname=self.filename 178 else: 179 self.nickname=nickname 180 f=pyfits.open(self.filename) 181 self.f=f 182 h=f[0].header 183 self.h=h 184 try: #Assume data is in extension 1. if not, fall back to extension 0. 185 self.data=f[1].data #.astype('Float32') 186 self.extnum=1 187 except IndexError: 188 self.data=f[0].data 189 self.extnum=0 190 self.exptime=h['exptime'] 191 self.camera=h['camera'] 192 self.saa_time=h['saa_time'] 193 self.badfile=osfn(h['maskfile']) 194 self.tdkfile=osfn(h.get('saadfile')) 195 self.gainplot=h['adcgain'] 196 197 self.inq1=slice(10,118),slice(10,118) 198 self.inq2=slice(10+128,118+128),slice(10,118) 199 self.inq3=slice(10+128,118+128),slice(10+128,118+128) 200 self.inq4=slice(10,118),slice(10+128,118+128) 201 202 self.q1=slice(0,128),slice(0,128) 203 self.q2=slice(128,256),slice(0,128) 204 self.q3=slice(128,256),slice(128,256) 205 self.q4=slice(0,128),slice(128,256) 206 207 208 print self.nickname, ": using DQ extension for badpix" 209 try: 210 self.dq=f['dq',1].data 211 if self.dq is not None: 212 dqmask = 1+16+32+64+128+256 #selected values 213 self.nonsourcemask=N.bitwise_and(self.dq, dqmask+1024) #exclude sources 214 self.nonsourceidx=N.where(self.nonsourcemask == 0) 215 self.nonsource=self.data[self.nonsourceidx] 216 self.badpix=N.bitwise_and(self.dq, dqmask) 217 else: 218 self.badpix=None 219 except KeyError,e: 220 print e 221 print 'DQ extension not found for %s'%imgfile 222 print 'defaulting to maskfile' 223 self.badpix = None 224 225 # print "self.badpix.shape = ",self.badpix.shape 226 227 if self.badpix is None: 228 print "failing over to ",self.badfile 229 try: 230 f2=pyfits.open(self.badfile) 231 self.badpix=f2['dq',1].data 232 f2.close() 233 except IOError,e: 234 print e 235 print "Bad pixel image not read" 236 print "Bad pixel image filename obtained from ",self.filename 237 self.badpix=None 238 239 #Don't leave file handles hanging around 240 self.f.close()

241

242 - def writeto(self,outname,clobber=False):

243 f=pyfits.open(self.filename) 244 f[self.extnum].data=self.data 245 f[0].header=self.h #update the primary header 246 f.writeto(outname,clobber=clobber)

247

248 - def dark_subtract(self,dark):

249 self.data=(self.data-dark)/self.exptime

250 251

252 - def pedskyish(self):

253 """ Performs something like the IRAF pedsky task, but with a bit more 254 sophistication in handling the central row and column""" 255 256 #Compute the median for each quadrant independently 257 ## m=N.array([imstat(self.data[self.inq1],nclip=0,binwidth=0.01,fields='median').median, 258 ## imstat(self.data[self.inq2],nclip=0,binwidth=0.01,fields='median').median, 259 ## imstat(self.data[self.inq3],nclip=0,binwidth=0.01,fields='median').median, 260 ## imstat(self.data[self.inq4],nclip=0,binwidth=0.01,fields='median').median]) 261 262 m=N.array([median(self.data[self.inq1]), 263 median(self.data[self.inq2]), 264 median(self.data[self.inq3]), 265 median(self.data[self.inq4])]) 266 ## print "file ",self.filename 267 ## print "raw m",m 268 ## temp=imstat(m,nclip=1,binwidth=0.01,fields='mean,median') 269 ## print "stats: mean/median",m.mean(),median(m) 270 m=m-median(m) 271 ## print "after sub",m 272 273 #Subtract the median from each quadrant 274 self.data[self.q1]=self.data[self.q1]-m[0] 275 self.data[self.q2]=self.data[self.q2]-m[1] 276 self.data[self.q3]=self.data[self.q3]-m[2] 277 self.data[self.q4]=self.data[self.q4]-m[3] 278 279 280 #"special handling of middle col/row" 281 if self.camera < 3: 282 temp=imstat( (self.data[:,127]-self.data[:,126]), 283 nclip=1,binwidth=0.01,fields='midpt') 284 self.data[:,127]=self.data[:,127]-temp.midpt 285 elif self.camera==3: 286 temp=imstat( (self.data[127,:]-self.data[126,:]), 287 nclip=1,binwidth=0.01,fields='midpt') 288 self.data[127,:]=self.data[127,:]-temp.midpt 289 else: 290 raise ValueError, "Bad camera value"

291 292 ##................................................................................... 293 ## Original code that I think is wrong: 294 ## transcribed parens from idl code incorrectly 295 ##................................................................................... 296 ## #Camera 3 is special: treat its middle column in a similar way 297 ## if self.camera < 3: 298 ## temp=imstat(self.data[:,127],nclip=1,binwidth=0.01,fields='median') 299 ## # print "line 127 median is ",temp.median 300 ## # print "line 127 mean is ",self.data[:,127].mean() 301 ## self.data[:,127]=self.data[:,127]-temp.median-self.data[:,126] 302 ## elif self.camera==3: 303 ## temp=imstat(self.data[127,:],nclip=1,binwidth=0.01,fields='median') 304 ## self.data[127,:]=self.data[127,:]-temp.median-self.data[126,:] 305 ## else: 306 ## raise ValueError, "Bad camera value" 307 308

309 - def getmask(self,dim=256,border=3,writename='mask.dat',clobber=False):

310 """Computes a mask to use for pixels to omit""" 311 mask=N.zeros((dim,dim),dtype=N.dtype('float32')) 312 badmask=N.ones((dim,dim),dtype=N.dtype('float32')) 313 if self.badpix is not None: 314 u=N.where(self.badpix != 0) 315 mask[u]=1 316 badmask[u]=0 317 # Always Mask out central "cross" chipgap 318 mask[(dim/2)-1,:]=1 319 mask[:,(dim/2)-1]=1 320 # and the very edges 321 mask[0:16,:]=1 #apparently the bottom edge is different 322 mask[dim-border:dim,:]=1 323 mask[:,0:border+1]=1 324 mask[:,dim-border:dim]=1 325 326 if writename: 327 writeimage(mask,writename,clobber=clobber) 328 return mask,badmask

329

330 - def apply_mask(self,mask):

331 goodpix=N.where(mask == 0) 332 self.masked_data = self.data[goodpix]

333 334

335 - def getscales(self,saaper,mask,pars):

336 337 cal=self.data*self.exptime 338 acc=saaper*self.exptime 339 340 for dom in self.domains.values(): 341 try: 342 sz1=int(dom.range/pars.stepsize)+1 343 stepval=[pars.stepsize*i for i in xrange(sz1)] 344 345 #there's got to be a better way to do this! 346 #Make a mask & fill it all with ones 347 fitmask=N.ones(mask.shape) 348 #Then make the pixels we want be set to zero 349 fitmask[dom.pixlist]=0 350 #Then set the mask-defined bad pixels to one so we don't use them 351 #(Notice there's no use of "self.badpix" here, wonder why not?) 352 #Ah! It's because self.badpix was already used in *making* that mask. OK. 353 badpix=N.where(mask == 1) 354 fitmask[badpix]=1 355 #Finally, choose only those pixels where it's set to zero. 356 umask=N.where(fitmask == 0) 357 358 dom.pp=N.zeros((3,int(dom.range/pars.stepsize)+1),dtype=N.dtype('float32')) 359 index=0 360 for i in stepval: 361 dif=cal-(acc*i) 362 temp=imstat(dif[umask],binwidth=0.01,nclip=3,fields='stddev,mode') #sigma=100 363 dom.pp[:,index]=i,temp.stddev,temp.mode 364 index+=1 365 dom.striplowerthan(pars.binsigfrac) 366 if pars.diagfile: 367 dom.writeto(pars.diagfile+'_'+dom.name+'_signal_domain.dat',clobber=pars.clobber) 368 ubest,umode=dom.getmin() 369 best=dom.pp[0,ubest] 370 371 print "\nResults summary for %s domain:"%dom.name 372 if pars.dofit: 373 minx=max(ubest-5,0) 374 maxx=min(ubest+5,len(dom.pp[0])-1) 375 thedata=[(dom.pp[0,i],dom.pp[1,i]) for i in range(minx,maxx+1)] 376 377 best,dom.chi2,itertrace=parabola_min(thedata,best) 378 # best=parabola1(dom.pp[0,minx:maxx],pp[1,minx:maxx],minguess=best) 379 # best=parabola1(dom.pp[0,minx:maxx],pp[1,minx:maxx],minguess=best) 380 381 dom.nr=(1.0-dom.pp[1,ubest]/dom.pp[1,0])*100 382 dom.scale=best 383 dom.bestloc=ubest 384 385 386 #print " zero-mode scale factor is : ",dom.pp[0,umode] 387 print " min-noise (best) scale factor is: ",dom.scale 388 print " effective noise at this factor (electrons at gain %f): %f"%(self.gainplot,dom.pp[1,ubest]*self.gainplot) 389 print " noise reduction (percent) : ",dom.nr 390 391 #Apply a sensibility check 392 if dom.scale < 0: 393 raise NegScaleError, "ERROR: Best scale factor for %s domain is negative"%dom.name 394 except ValueError,e: 395 print "Error calculating scale for %s domain"%dom.name 396 print str(e) 397 print "No correction can be calculated for this domain" 398 dom.nr = 0 399 dom.scale = 0 400 dom.bestloc=0 401 dom.chi2=0

402

403 - def apply_domains(self,saaper,badmask,noisethresh,appimage=None):

404 if appimage is not None: 405 final=appimage 406 else: 407 final=self.data.copy() 408 409 saacorr=N.zeros(final.shape,dtype=N.dtype('float32')) 410 411 hdom,ldom=self.domains['high'],self.domains['low'] 412 self.update=1 413 if hdom.nr >= noisethresh and ldom.nr >= noisethresh: 414 print "\n Applying noise reduction in both domains " 415 self.appstring='both' 416 saacorr[ldom.pixlist]=saaper[ldom.pixlist]*(ldom.scale*badmask[ldom.pixlist]) 417 saacorr[hdom.pixlist]=saaper[hdom.pixlist]*(hdom.scale*badmask[hdom.pixlist]) 418 419 elif hdom.nr > noisethresh and ldom.nr < noisethresh: 420 print "\n Applying noise reduction in high domain only " 421 self.appstring='high only' 422 saacorr[hdom.pixlist]=saaper[hdom.pixlist]*(hdom.scale*badmask[hdom.pixlist]) 423 424 elif hdom.nr < noisethresh and ldom.nr >= noisethresh: 425 print "\n...Noise reduction in high domain < 1%: applying low scale everywhere" 426 self.appstring='low everywhere' 427 saacorr=saaper*(ldom.scale*badmask) 428 429 elif hdom.nr < noisethresh and ldom.nr < noisethresh: 430 print "\n*** Noise reduction < 1 %, not applying" 431 self.appstring='none' 432 self.update=0 433 else: 434 raise ValueError,"Huh?? hi_nr, lo_nr: %f %f"%(hdom.nr,ldom.nr) 435 436 if self.appstring != 'none': 437 final=final-saacorr 438 439 ## import futil 440 ## futil.writeimage(saacorr,'scaled_sapper.fits') 441 442 return final

443

444 - def update_header(self,pars,tag='default',header=None):

445 """ Update the FITS header with all this good stuff we've done""" 446 447 #Start with the last keyword, for ease of applying. 448 449 if header is None: 450 header=self.h 451 452 #Describe what was applied 453 lastkey='SCNAPPLD' 454 header.update(lastkey, 455 self.appstring, 456 'to which domains was SAA cleaning applied', 457 after='SAACRMAP') 458 459 460 #Then work forward from the beginning of the section: 461 462 #First put in a comment card as a separator 463 header.add_blank('',before=lastkey) 464 header.add_blank(' / SAA_CLEAN output keywords',before=lastkey) 465 header.add_blank('',before=lastkey) 466 467 #Then describe the persistence image: 468 header.update('SAAPERS', 469 os.path.basename(pars.saaperfile), 470 'SAA persistence image', 471 before=lastkey) 472 if not pars.readsaaper: 473 header.update('SCNPSCL', 474 pars.scale, 475 'scale factor used to construct persistence img', 476 before=lastkey) 477 header.update('SCNPMDN', 478 pars.saaper_median, 479 'median used in flatfielding persistence image', 480 before=lastkey) 481 header.add_blank('',before=lastkey) 482 483 #Describe the domains 484 header.update('SCNTHRSH', 485 self.thresh, 486 'Threshold dividing high & low signal domains', 487 before=lastkey) 488 header.update('SCNHNPIX', 489 self.domains['high'].npix, 490 'Number of pixels in high signal domain (HSD)', 491 before=lastkey) 492 header.update('SCNLNPIX', 493 self.domains['low'].npix, 494 'Number of pixels in low signal domain (LSD)', 495 before=lastkey) 496 header.add_blank('',before=lastkey) 497 498 #Describe the results in each domain 499 ## self.h.update('SCNGAIN', 500 ## self.gainplot, 501 ## 'gain used for effective noise calculations', 502 ## before=lastkey) 503 for k in self.domains: 504 HorL=k[0].upper() 505 header.update('SCN%sCHI2'%HorL, 506 self.domains[k].chi2, 507 '%sSD chi squared for parabola fit'%HorL, 508 before=lastkey) 509 header.update('SCN%sSCL'%HorL, 510 self.domains[k].scale, 511 '%sSD scale factor for min noise'%HorL, 512 before=lastkey) 513 bestloc=self.domains[k].bestloc 514 header.update('SCN%sEFFN'%HorL, 515 self.domains[k].pp[1,bestloc]*self.gainplot, 516 '%sSD effective noise at SCNGAIN'%HorL, 517 before=lastkey) 518 header.update('SCN%sNRED'%HorL, 519 self.domains[k].nr, 520 '%sSD noise reduction (percent)'%HorL, 521 before=lastkey)

522 ##................................................................ 523 ## Only needed for testing: removed for release 524 ##................................................................ 525 ## self.h.update('SCNTAG', 526 ## tag, 527 ## 'Tag/description of this version', 528 ## 529 ##................................................................ 530 531 #.......................................................................... 532 # Exception definitions

533 -class NoPersistError(exceptions.Exception):

534 pass

535 -class BadThreshError(exceptions.Exception):

536 pass

537 -class NegScaleError(exceptions.Exception):

538 pass

539 -class InsuffImprovement(exceptions.Exception):

540 pass

541 -class AlreadyDone(exceptions.Exception):

542 pass

543 #............................................................................. 544 #Helper functions: 545 #-............................................................................

546 -def osfn(filename):

547 """Return a filename with iraf syntax and os environment names substituted out""" 548 if filename is None: 549 return filename 550 551 #Baby assumptions: suppose that the env variables will be in front. 552 553 if filename.startswith('$'): #we need to translate a logical 554 symbol,rest=filename.split('/',1) 555 elif '$' in filename: #we need to fix the iraf syntax 556 symbol,rest=filename.split('$',1) 557 else: 558 return filename 559 newfilename=os.environ[symbol]+'/'+rest 560 return newfilename

561

562 -def writeimage(image, filename, comment=None,clobber=False):

563 hdulist=pyfits.HDUList() 564 hdu=pyfits.PrimaryHDU() 565 hdu.data=image 566 if (comment is not None): 567 hdu.header.add_comment(comment) 568 hdulist.append(hdu) 569 hdulist.writeto(filename,clobber=clobber)

570 571 #.......................................................................... 572 # Math functions

573 -def median(a):

574 return N.sort(a.ravel())[a.size / 2]

575

576 -def parabola_model(coeffs,t):

577 r=coeffs[0]*(t-coeffs[1])**2 + coeffs[2] 578 return r

579

580 -def parabola_min(thedata, startguess):

581 #We may not need to rescale the data 582 guesscoeff=(100,startguess,0.1) 583 fitcoeff,chi2,itertrace=LeastSquares.leastSquaresFit(parabola_model,guesscoeff,thedata) 584 print "chi2 for parabola fit = ",chi2 585 return fitcoeff[1],chi2,itertrace

586 #............................................................................ 587 #Functions for gauss-poly fitting of the persistence image histogram

588 -def gausspoly_eval(coeffs,t):

589 z=(t-coeffs[1])/coeffs[2] 590 zz=-1*(z**2/2.) 591 r=coeffs[0]*N.exp(zz) + coeffs[3] + coeffs[4]*t + coeffs[5]*t**2 592 return r

593

594 -def gausspoly_model(coeffs,t):

595 import math 596 z=(t-coeffs[1])/coeffs[2] 597 zz=-1*(z**2/2.) 598 r=coeffs[0]*zz.exp() + coeffs[3] + coeffs[4]*t + coeffs[5]*t**2 599 return r

600

601 -def gausspoly_fit(thedata,guesscoeff):

602 #coeffs are: 603 # amplitude of gaussian: guess max(data) 604 # center of gaussian: guess mode(data) 605 # sigma of gaussian: guess stddev(data) 606 # constant: guess 0.1 607 # linear term: guess 0.1 608 # quadratic term: guess 0.0 609 610 fitcoeff,chi2,itertrace=LeastSquares.leastSquaresFit(gausspoly_model, 611 guesscoeff, 612 thedata) 613 return fitcoeff, chi2, itertrace

614 615 #...........................................................................

616 -def thresh_from_gausspoly_fit(saa, parbinwidth=0.5, nclip=3, 617 diagfile=None, clobber=False):

618 """ Some massaging of the SAApersistence image histogram is 619 performed in order to obtain an optimal fit. 620 Unfortunately this involves some magic numbers taken from 621 the IDL code.""" 622 #Work in a bigger dynamic range space 623 im=saa*500.0 624 binwidth=parbinwidth*500.0 625 626 #Compute the histogram 627 hnbins=int( (10000+100)/binwidth) + 1 628 h=histogram1d(im,hnbins,binwidth,-100) 629 xloc=N.arange(h.nbins)*h.binWidth+h.minValue 630 631 #Select out only the data range we're interested in 632 #Take a first guess at the standard deviation 633 idx=((im >= -100) & (im <= 10000)) 634 yy=imstat(im[idx],binwidth=0.1,nclip=3,fields='stddev') 635 636 #Set up the data we're going to fit 637 if hnbins > 600: 638 numpoints=600 639 else: 640 numpoints=hnbins 641 642 thedata = [(xloc[i],h.histogram[i]) for i in range(numpoints)] 643 t=xloc[0:numpoints] 644 645 #Now set up the start guesses 646 hmax=h.histogram[0:numpoints].max() 647 hbinmax=xloc[h.histogram[0:numpoints].argmax()] 648 startguess=[hmax, 649 xloc[h.histogram[0:numpoints].argmax()], 650 yy.stddev, 651 0.1, 0.1, 0.0] 652 653 #Do the fitting: with a catch for a linear algebra failure 654 try: 655 coeffs,chi2,itertrace=gausspoly_fit(thedata,startguess) 656 except LinAlgError, e: 657 if diagfile is None: 658 diagfile='diag_linalgerr' 659 f=smartopen(diagfile+'_gp_hist.txt','w',clobber=clobber) 660 for k in range(len(t)): 661 line = "%f %d\n"%(t[k],h.histogram[k]) 662 f.write(line) 663 f.close() 664 raise e 665 666 #and tell us about the results 667 print "\nCoefficients for gauss-poly fit to persistence model histogram:" 668 r=itertrace[-1] #Last iteration 669 670 print "Gaussian (low signal component) terms:" 671 print " Amplitude, Mean, Sigma: %f %f %f"%(r[0].value,r[1].value,r[2].value) 672 print "Polynomial terms:" 673 print " Constant, Linear, Quadratic:%f %f %f"%(r[3].value,r[4].value,r[5].value) 674 print"" 675 676 if diagfile: 677 678 #This prints the histogram that is actually fit 679 f=smartopen(diagfile+'_gp_hist.txt','w',clobber=clobber) 680 for k in range(len(t)): 681 line = "%f %d\n"%(t[k],h.histogram[k]) 682 f.write(line) 683 f.close() 684 685 #and this prints the series of coefficients for each iteration 686 f=smartopen(diagfile+'_gp_iters.txt','w',clobber=clobber) 687 for p in (itertrace): 688 line=' '.join([str(x[0]) for x in p])+"\n" 689 f.write(line) 690 f.close() 691 692 ## for k in range(len(itertrace)): 693 ## print "Iter %d"%k,itertrace[k] 694 695 #Finally, compute the threshold based on the fit. 696 #Don't forget to divide out the magic-number to convert back 697 #to the original scale. 698 699 thresh=(coeffs[1] + 3.5*abs(coeffs[2]))/500. 700 return thresh

701 702 #.............................................................................. 703 # General functions 704 #..........................................................................

705 -def get_postsaa_darks(imgfile):

706 """ Return the filenames containing the post-saa dark exposures, if 707 present. Otherwise raise an exception and exit. """ 708 709 #Get the science header 710 inpath=os.path.dirname(osfn(imgfile)) 711 if inpath != '': 712 inpath+= '/' 713 f=pyfits.open(imgfile) 714 h=f[0].header 715 saa_asn=h['saa_dark'] 716 f.close() 717 if saa_asn == 'N/A': 718 raise NoPersistError, """This data was not taken in an SAA-impacted orbit. 719 No correction needed. Exiting.""" 720 else: 721 #Get the files out of that set 722 saa_files=[] 723 f2=pyfits.open(inpath+saa_asn.lower()+'_asn.fits') 724 for i in [0,1]: 725 name=f2[1].data[i] 726 saa_files.append(inpath+name.field(0).lower()+'_raw.fits') 727 f2.close() 728 return saa_files

729

730 -def getdark(camera,tdkfile,darkpath):

731 """ Get the right dark file for a given NICMOS camera. 732 This is definitely not the right way to do this.""" 733 dfile={1:'c1_saadarkref_drk.fits', 734 2:'c2_saadarkref_drk.fits', 735 3:'c3_saadarkref_drk.fits'} 736 darkpath=os.path.abspath(darkpath)+'/' 737 defaultfile=darkpath+dfile[camera] 738 #Choose which file to use 739 if tdkfile: 740 altfile=darkpath+os.path.basename(tdkfile) 741 if os.path.isfile(tdkfile): 742 thefile=tdkfile 743 elif os.path.isfile(altfile): 744 thefile=altfile 745 else: 746 thefile=defaultfile 747 else: 748 thefile=defaultfile 749 #Return the data from it 750 f=pyfits.open(thefile) 751 ans= f[1].data 752 f.close() 753 return ans

754

755 -def make_saaper(im1,im2,dark,pars,crthresh=1):

756 #Process the data 757 for im in [im1,im2]: 758 im.dark_subtract(dark) 759 im.pedskyish() 760 #Combine the data 761 saaper=((im1.data*pars.wf1) + (im2.data/pars.scale)*pars.wf2) 762 #Correct for CRs 763 if pars.crthresh: 764 a=im1.data-(im2.data/pars.scale) 765 u1=N.where(a > pars.crthresh) 766 saaper[u1]=im2.data[u1]/pars.scale 767 768 a=(im2.data/pars.scale) - im1.data 769 u2=N.where(a > pars.crthresh) 770 saaper[u2]=im1.data[u2] 771 if pars.writesaaper and pars.saaperfile: 772 writeimage(saaper,pars.saaperfile,clobber=pars.clobber) 773 return saaper

774 775

776 -def get_dark_data(imgfile,darkpath):

777 saafiles=get_postsaa_darks(imgfile) 778 im1=Exposure(saafiles[0],nickname='postsaa dark #1') 779 im2=Exposure(saafiles[1],nickname='postsaa dark #2') 780 dark=getdark(im1.camera,im1.tdkfile,darkpath) 781 return im1,im2,dark

782

783 -def flat_saaper(saaper,img):

784 #The "midpt" computes a pseudo-median function. 785 mm=imstat(saaper,nclip=1,binwidth=0.01,fields='midpt').midpt 786 #Use median, or mode? which is better? 787 if img.h['flatdone'] == 'PERFORMED': 788 flatname=osfn(img.h['flatfile']) 789 ## if flatname.startswith('nref$'): 790 ## prefix,root=flatname.split('$',1) 791 ## flatname=iraf.osfn(prefix+'$')+root 792 flat=Exposure(flatname,nickname='flatfile') 793 794 print "median used in flatfielding: ",mm 795 saaper=((saaper-mm)*flat.data) + mm 796 return saaper,mm

797 798

799 -def smartopen(fname, mode, clobber=True):

800 """ Allows specifying a clobber behavior """ 801 if mode.startswith('w') and not clobber: 802 if os.path.isfile(fname): 803 raise IOError, "%s already exists"%fname 804 805 handle=open(fname,mode) 806 return handle

807 808 #.................................................................... 809 # The "main" program 810 #....................................................................

811 -def clean(usr_calcfile,usr_targfile,usr_outfile,pars=None):

812 numerixenv.check() #Temporary NUMERIX environment check 813 print "Input files: %s %s"%(usr_calcfile,usr_targfile) 814 imgfile=osfn(usr_calcfile) 815 img=Exposure(imgfile,nickname='sci image') 816 targfile=osfn(usr_targfile) 817 if imgfile != targfile: 818 #then we'll need the data from the targfile 819 targ=Exposure(targfile,nickname='target image') 820 appimage=targ.data.copy() 821 else: 822 #we'll apply it to the same file 823 targ=img 824 appimage=None 825 826 827 #Check to make sure we're not trying to run the task twice 828 #on the same image 829 already_done=['low everywhere','both','high only'] 830 for check in [targ,img]: 831 scnappld=check.h.get('scnappld',None) 832 if scnappld in already_done: 833 raise AlreadyDone, check.filename 834 835 outfile=osfn(usr_outfile) 836 if pars is None: 837 pars=params() 838 if pars.readsaaper: 839 sfile=pyfits.open(pars.saaperfile) 840 saaper=sfile[0].data 841 sfile.close() 842 else: 843 im1,im2,dark=get_dark_data(imgfile,pars.darkpath) 844 saaper=make_saaper(im1,im2,dark,pars) 845 print "Using scale factor of ",pars.scale," to construct persistence image" 846 847 848 mask,badmask=img.getmask(writename=pars.maskfile,clobber=pars.clobber) 849 saaper,mm=flat_saaper(saaper,img) 850 pars.saaper_median=mm 851 852 if pars.flatsaaperfile: 853 writeimage(saaper,pars.flatsaaperfile,clobber=pars.clobber) 854 855 mask,badmask=img.getmask(writename=None) 856 img.apply_mask(mask) 857 858 if pars.fitthresh: 859 img.thresh = thresh_from_gausspoly_fit(saaper, 860 parbinwidth=pars.histbinwidth, 861 nclip=pars.nclip, 862 diagfile=pars.diagfile, 863 clobber=pars.clobber) 864 else: 865 img.thresh=pars.thresh 866 867 868 print "Threshold for hi/lo: ",img.thresh 869 870 #Apply threshold *to persistence image* 871 img.domains={'high':Domain('high', 872 N.where(saaper > img.thresh), 873 pars.hirange), 874 'low' :Domain('low', 875 N.where(saaper <= img.thresh), 876 pars.lorange) 877 } 878 879 880 print "Npixels hi/lo: ",len(img.domains['high'].pixlist[0]),len(img.domains['low'].pixlist[0]) 881 882 #Do some checking for sensible results 883 ## if (img.domains['high'].npix == 0): 884 ## raise BadThreshError,"ERROR: Zero pixels found in high signal domain." 885 if (img.domains['high'].npix > img.domains['low'].npix): 886 raise BadThreshError,"ERROR: Number of high domain pixels exceeds the number of low domain pixels" 887 img.getscales(saaper,mask,pars) 888 889 final=img.apply_domains(saaper,badmask,pars.noisethresh,appimage=appimage) 890 891 892 if 1: #img.update: 893 targ.data=final 894 img.update_header(pars,header=targ.h) 895 targ.writeto(outfile,clobber=pars.clobber) 896 897 898 return saaper,img

899

Source Code for Module nictools.saaclean