;+ ; NAME: ; make_spe_zero1 ; ; PURPOSE: (one line) ; make zero-flux (dark+bias) averaged files from spe files ; ; DESCRIPTION: ; Given a filenames of SPE zero-flux files, make averaged zero-flux file ; ; CATEGORY: ; IMAGES ; ; CALLING SEQUENCE: ; make_spe_zero1,filename, avg, err, bad, d = d, $ ; silent=silent, median = median, $ ; imbad = imbad, gain=gain, bias=bias, minbias=minbias, $ ; readn_dn = readn_dn, $ ; caldir=caldir, outfilename = outfilename, $ ; overrideoutfilename = overrideoutfilename ; INPUTS: ; filename : filename of SPE image cube with zero-flux (dark+bias) images ; ; OPTIONAL INPUT PARAMETERS: ; -------- behavior of the routine ------------ ; silent : 1 if the routine is to perform silently. Default = 0, ; print out statements ; median : 1 if using median averaging, 0 if using avgclip algorithm ; ; -------- description of the input file ------ ; imbad : if set, an array of bad image indices (eg ; imbad=[0] to drop 1st image ; bias, gain, and readn_dn are used in some modes to calculate the ; expected noise per pixel, used for flagging outliers ; bias : if set, an estimate of the bias level. Default = 0 or the ; minimum of the input file (excluding bad images) ; minbias : if set, default bias is the minimum of the input file ; gain : electron/dn. Default = 1. ; readn_dn : read noise in dn. Default = caculate from the data ; exptime : if set, the exposure time in seconds, otherwise ; exposure time is taken from the SPE header (needed for ; triggered zero-flux files with PhotonMax systems) ; ; -------- output file ------------ ; caldir : data directory for cal files (with terminal '/'). Default = 'cal/' ; overrideoutfilename : if set (1), use outfilename as output file name. ; Default = 0, construct the output filename from the input file name ; outputfilename : if overrideoutfilename is set, this is the output ; filename used ; exptime : if set, use passed value for EXPTIME keyword ; ; OUTPUTS: ; avg - the average zero-flux image ; err - the error in the average zero-flux image (NOT error per ; image per pixel) ; bad - bad pixel image for avg (0 = good) ; ; OPTIONAL OUTPUTS: ; d - the input file (after deselecting with imbad) ; ; SIDE EFFECTS: ; For each file, writes ; cal/_.fits ; ; For each combination of gain, adc, and binning, write ; /_zero.fits ; with ; primary HDU - ; - INFILE - input file used to make average ; - OBJECT - 'median image of zero-flux files' ; - NIMAGES - number of images used to create average ; - EXPTIME - Exposure time in seconds ; - READN_DN - Derived noise per pixel, read noise + dark noise (DN) ; - AVG_METH - Method for making the average ; - header converted from original SPE header ; extention 1 - zero-flux error DN ; extention 2 - bad pixel map (0=good) ; ; EXAMPLE: ; ; MODIFICATION HISTORY: ; Written 14 Jun 2007 Leslie Ann Young SwRI ; Modified 13 Sep 2007 LAY ; Modified 18 Sep 2007 LAY added exptime keyword ; Modified 20 Sep 2007 LAY print input file ; Modified 21 Sep 2007 LAY median only the middle 100 frames ;- pro make_spe_zero1,filename, avg, err, bad, d = d, $ silent=silent, median = median, $ imbad = imbad, gain=gain, bias=bias, minbias=minbias, $ readn_dn = readn_dn, $ caldir=caldir, outfilename = outfilename, $ overrideoutfilename = overrideoutfilename, $ exptime = exptime if not keyword_set(silent) then silent = 0b ;-------------------------------- ; Read the file using filename ;-------------------------------- fullfn = filename[0] if not silent then begin print, 'make_spe_zero1: averaging ', fullfn endif fn = file_basename(fullfn) len = strlen(fn) if strmatch(fn,'*.gz', /fold) then begin basefn = strmid(fn,0,len-7) endif else begin basefn = strmid(fn,0,len-4) endelse read_princeton_gz,fullfn,d,head=h ; n_elements rather than keyword_set in case ; user passes exptime = 0 if n_elements(exptime) eq 0 then begin exptime = h.exp_sec endif if not silent then begin print, 'EXPTIME = ', float(exptime) end nxnyn, d, nx, ny, nimages ;-------------------------------- ; Filter out the bad images ;-------------------------------- if n_elements(imbad) gt 0 then begin if not silent then begin print, 'DESELECTING IMAGES ', imbad end isgood = bytarr(nimages) + 1b isgood[imbad] = 0 igood = where(isgood eq 1b, ngood) if ngood eq 0 then begin print, 'make_spe_zero1(', filename,'...): all images deselected' return endif else begin d = d[*,*, igood] endelse nxnyn, d, nx, ny, nimages end ;-------------------------------- ; Set readnoise (rn). ; If the readnoise is not specified, estimate it from the ; difference in the first two images ;-------------------------------- if keyword_set(readn_dn) then begin rn = readn_dn endif else begin if nimages eq 1 then begin d0 = float(d[*,*,0]) robomean,d0,5.,0.1,avg,avgdev,stddev,var,skew,kurt,nfinal rn = stddev endif else begin d0 = float(d[*,*,1]) - float(d[*,*,0]) robomean,d0,5.,0.1,avg,avgdev,stddev,var,skew,kurt,nfinal rn = stddev / sqrt(2.) endelse endelse if not silent then begin print, 'READNOISE (1st pass) = ', rn, 'DN' end ;-------------------------------- ; Set bias (b). ; If the bias is not specified, set to zero or ; the minimum of the input file ;-------------------------------- if keyword_set(bias) then begin b = bias endif else begin if keyword_set(minbias) then begin b = min(d) endif else begin b = 0. endelse endelse if not silent then begin print, 'MINVAL = ', b, 'DN' end ;-------------------------------- ; Set gain (g). ; If the gain is not specified, set to 1 ;-------------------------------- if keyword_set(gain) then begin g = gain endif else begin g = 1. endelse if not silent then begin print, 'ASSUMED GAIN = ', g, 'e/DN' end ;-------------------------------- ; Get the zero/bias frame ; ;-------------------------------- if keyword_set(median) then method = 'median' else method = 'avgclipsim' case nimages of ;-------------------------------- ; 1 image: ; avg = that image ; err = passed read noise, or estimated read noise from robomean ; bad = 1 where points are more than 5 sigma from a 10x10 median-smoothed ; version of the average ;-------------------------------- 1: begin thresh = 5.0 avg = d[*,*,0] err = replicate(rn, nx, ny) ; error in the zero smavg = median(avg, 10) bad = bytarr(nx,ny) indxbad = where( abs(avg - smavg)/err gt thresh, nbad) if nbad gt 0 then bad[indxbad] = 1b end ;-------------------------------- ; 2 images: ; avg = mean of the two images ; err = passed read noise, or estimated read noise from robomean ; bad = 1 where the difference is > 5 sigma (5 * rn * sqrt(2)) ;-------------------------------- 2: begin thresh = 5. avg = (d[*,*,0] + d[*,*,1])/2. err = replicate(rn, nx, ny)/2. ; error in the zero bad = bytarr(nx,ny) indxbad = where (abs(d[*,*,0] - d[*,*,1])/(sqrt(2.)*rn) gt thresh, nbad) if nbad gt 0 then bad[indxbad] = 1b end else: begin case method of ;-------------------------------- ; >=3 images, median: ; avg = median ; err = std deviation of mean sqrt( total(resid^2)/( (n-1)*n) ) ; bad = 1 where error is 5 sigma larger than other errors ;-------------------------------- 'median': begin ; call the (slow) medarr_mwb rather ; than risk malloc error if n_images gt 1000 then begin medarr_mwb, d, avg endif else begin avg = median(d, dim=3, /even) endelse ; create the error array res = d for i = 0, nimages-1 do res[*,*,i] = d[*,*,i]-avg ; err1 is the error in one pixel (i.e., the readnoise) err1 = fltarr(nx, ny) for i=0, nx-1 do begin for j = 0, ny-1 do begin err1[i,j] = total(res[i,j,*]^2) endfor endfor err1 = sqrt( err1/(nimages-1) ) ; error in a single pixel err = err1/sqrt(nimages) ; error in the average robomean,err1,5.,0.1,avgerr,avgerrdev,stddeverr print, 'READNOISE = ', avgerr, 'DN' ; create the bad-pixel array bad = bytarr(nx, ny) ibad = where(err gt (avgerr + 5.0 * stddeverr), nbad) if nbad gt 0 then bad[ibad] = 1 end ; median ;-------------------------------- ; >=3 images, median: ; avg = average, clipping 3-sigma outliers ; err = std deviation of mean excluding marked outliers ; bad = 1 where all pixels are outliers ;-------------------------------- 'avgclipsim': begin thresh = 3.0 ; 2007 Sept 21: malloc errors when ; median'ing 0.5-Gb file if n_images gt 1000 then begin first = (nimages / 2 - 500) > 0 last = (first + 999) < (nimages - 1) avg = median(d[first:last], dim=3, /even) endif else begin avg = median(d, dim=3, /even) endelse ; approximate noise, avgclip algorithm sigma = ( sqrt( ( (avg-b) > 0.)/g ) > rn ) badmat = bytarr(nx,ny,nimages) asum = fltarr(nx,ny) acnt = fltarr(nx,ny) npts = long(nx) * long(ny) npixclipped = 0 for i=0,nimages-1 do begin new = d[*,*,i] resid = (new-avg)/sigma z = where(abs(resid) ge thresh,count) if count ne 0 then begin badi = badmat[*,*,i] badi[z] = 1B badmat[*,*,i] = badi endif z = where(badmat[*,*,i] eq 0B,count) if count ne 0 then begin asum[z] = asum[z]+new[z] acnt[z] = acnt[z] + 1.0 endif else begin print,' Serious Warning! all pixels in image ',strn(i), $ ' were considered bad.' ;stop endelse npixclipped = npixclipped + (npts - count) endfor ; refine the average where it's possible to z = where(acnt ne 0,count) if count ne 0 then avg[z]= asum[z]/acnt[z] ; create the bad pixel array bad = (acnt le 1) ; create the array with the error in avg err1 = replicate(rn, nx, ny) for i = 0, nx-1 do begin for j=0, ny-1 do begin n = acnt[i,j] if n ge 2 then begin err1[i,j] = $ sqrt( total( $ (d[i,j,*] - avg[i,j])^2 * (1.-badmat[i,j,*])) / $ (n-1.) ) endif endfor endfor err = err1 / sqrt(n) robomean,err1,5.,0.1,avgerr,avgerrdev,stddeverr if not keyword_set(silent) then begin print, 'READNOISE = ', avgerr, 'DN' endif end ; avgclipsim endcase end end ;-------------------------------- ; Get or create the output filename ;-------------------------------- if not keyword_set(caldir) then caldir = 'cal/' if keyword_set(overridefilename) then begin zerofn = outfilename endif else begin zerofn = caldir + basefn+'_zero.fits' outfilename = zerofn endelse if not silent then print, zerofn ;--------------------------------------------- ; write it out ;--------------------------------------------- mkhdr, header, avg, /EXTEND header = spe_head2fits(h, infits = header) sxaddpar, header, 'INFILE', fullfn, before='FIRSTSPE' sxaddpar, header, 'OBJECT', 'median image of zero-flux files', before='FIRSTSPE' sxaddpar, header, 'NIMAGES', nimages, 'Number of images in original file', $ before='FIRSTSPE' sxaddpar, header, 'EXPTIME', exptime, 'Exposure time in seconds', $ before='FIRSTSPE' sxaddpar, header, 'READN_DN', avgerr, $ 'Derived noise per pixel, read noise + dark noise (DN)', $ before='FIRSTSPE' sxaddpar, header, 'AVG_METH', method, 'Method for making the average', $ before='FIRSTSPE' WRITEFITS, zerofn, avg, header mkhdr, header, err, /EXTEND sxaddpar, header, 'OBJECT', 'error in median image' WRITEFITS, zerofn, err, /append mkhdr, header, bad, /EXTEND sxaddpar, header, 'OBJECT', 'bad pixel list (0=good)' WRITEFITS, zerofn, bad, /append end