;+ ; NAME: ; coaddgimg ; ; PURPOSE: ; Stack all guider frames during a science exposure so that we have a time-integrated PSF. Writes an output coadded image file. ; Optional Use for DOS: computes the coadded image and extracts the seeing and transparency to add into DOS FITS header. ; ; CALLING SEQUENCE: ; coaddgimg, plateid, mjd, exposure, outdir=outdir, guidedir=guidedir, dos=dos, hdr=hdr, filename=filename, drpver=drpver ; ; INPUTS: ; plateid - plate id ; mjd - mjd of observations ; exposure - exposure number of observations ; ; OPTIONAL INPUTS: ; outdir - set this to indicate the output directory to store file ; guidedir - set this to indicate the guider image file directory ; filename - DOS output filename of the exposure being reduced ; hdr - header of DOS output filename ; ; OPTIONAL OUTPUTS: ; psfpar - contains the output from the psf parameter estimation ; errorcode - contains the errorcode for different failure mode ; errorcode = 0 : success processing ; errorcode = 1 : No valid guider images found within the interval of the science exposure. ; ; KEYWORDS: ; dos - set this to skip writing the coadded image, and proceed straight to parameter estimation ; drpver - the version of MaNGA DRP being used ; ; OUTPUTS: ; none ; ; PROCEDURES USED: ; ml_getenv ; headfits ; mrdfits ; mwrfits ; sxpar ; sxaddpar ; ; REVISION HISTORY: ; Jan-10-2013 - Written by Renbin Yan ; Mar-24-2014 - Generalized for DOS and DRP integration (B. Cherinka) ; Mar-27-2014 - Changed hdr to guidhdr, and added hdr keyword for DOS (B. Cherinka) ; May-??-2014 - Made more robust for mising guider images for exposures ; May-??-2014 - Added in transparency calculation (Renbin) ; Jul-10-2014 - Fixed bug that breaks for local run without guider images (B. Cherinka) ; Sep-09-2014 - Add safety check that dra/ddec defined in case first ; frame missing (D. Law) ; Mar-11-2015 - Added output keyword to hold psfparam output ; ;- pro coaddgimg, plateid, mjd, exposure, outdir=outdir, guidedir=guidedir, dos=dos, hdr=hdr, filename=filename, drpver=drpver, psfpar=psfpar, errorcode=errorcode on_error,0 compile_opt idl2 compile_opt idl2, hidden common subimageblock, stamp,errstamp,nx common flatinfo, flatfile_inmem, flat,hdr_flat common darkinfo, darkfile_inmem, dark,hdr_dark errorcode = 0 ; Define variables and directories platestr = string(plateid,format='(i0.0)') mjdstr = string(mjd,format='(i0.0)') guidedir = n_elements(guidedir) eq 0 ? ml_getenv('GCAM_DATA')+mjdstr+'/' : guidedir drpver = ~keyword_set(drpver) ? mangadrp_version(/simple) : strtrim(drpver,2) outdir = ~keyword_set(outdir) ? djs_filepath('', root_dir=ml_getenv('MANGA_SPECTRO_REDUX'), subdir=[drpver,platestr,mjdstr]) : outdir ; Grab exposure files expodir = ml_getenv('MANGA_SPECTRO_DATA')+mjdstr+'/' expostr = string(exposure,format='(i8.8)') files = 'sdR-??-'+expostr+'.fit*' expofiles=file_search(expodir,files,count=ct) if ct eq 0 then if ~keyword_set(dos) then message,'No files found for exposure:'+expostr+' in path:'+expodir else begin splog, 'No file found for exposure:'+expostr+' in path:'+expodir splog, 'Cannot compute seeing and transparency from co-added guider image' return endelse ; Grab associated guider images file = expofiles[0] sci_hdr =headfits(file) ; Apply any header fixes ml_mghdrfix,file,sci_hdr,/silent tai_beg = long64(sxpar(sci_hdr,'TAI-BEG')) tai_end = tai_beg+fix(sxpar(sci_hdr,'EXPTIME')) guider1=sxpar(sci_hdr,'GUIDER1') guidern=sxpar(sci_hdr,'GUIDERN') if long(exposure) eq 183437 then guidern = 'proc-gimg-0246.gits.gz' seqstart = fix(strmid(guider1,10,4)) seqlast = fix(strmid(guidern,10,4))+2 nimg = seqlast-seqstart+1 filestr = mjdstr+'-'+expostr ; Loop over guider images firstframe = 1 ; Safety flag if first proc frame missing firstfail=0 ; false by default n_validframe = 0 for i=0,nimg-1 do begin seqno = string(seqstart+i,format='(i4.4)') origin = 'gimg-'+seqno+'.fits*' exist1 = file_search(guidedir+origin,count=c1) if c1 ne 1 then begin splog,'WARNING: missing guider file:,',guidedir+origin splog,' Skipping it in the coadd. But please find out why it is missing!!!' Wait, 3 continue endif guidhdr=headfits(guidedir+origin) expbegin = date2tai(sxpar(guidhdr,'DATE-OBS')) expend = expbegin +fix(sxpar(guidhdr,'EXPTIME')) if expend lt tai_beg or expbegin gt tai_end then begin splog,' Frame:'+seqno+' does not overlap within the science exposure interval.' continue endif proc ='proc-gimg-'+seqno+'.fits*' exist2 = file_search(guidedir+proc,count=c2) if strcompress(sxpar(guidhdr,'IMAGETYP'),/remove) eq 'dark' then begin ; added by Kai Zhang for exposure 242039 MJD 57860 plate 9778 splog,'WARNING: dark frame '+seqno+' is included, skip' proc='proc-gimg-'+string(seqstart,format='(i4.4)')+'.fits*' continue endif if c2 ne 1 then begin ; error = 1 splog,'WARNING: missing processed guider file:,',guidedir+proc splog,' Not a big deal. but I cannot check the consistency of the decenter keywords.' firstfail=1 dra=0. ddec=0. endif else begin proc_hdr=headfits(guidedir+proc) if firstframe then begin dra=sxpar(proc_hdr,'DCNRA') ddec=sxpar(proc_hdr,'DCNDEC') endif else begin ; DRL- add safety check that first frame didn't fail if (~firstfail) then begin if dra ne sxpar(proc_hdr,'DCNRA') or ddec ne sxpar(proc_hdr,'DCNDEC') then print,'WARNING: DCNRA or DCNDEC changed during this science exposure. !!!' endif endelse endelse img=mrdfits(guidedir+origin,0,/unsigned,/silent) imgsz = size(img,/dimen) if imgsz[0] ne 524 or imgsz[1] ne 512 then begin message,'Error: I am expecting the guider frames to be 524x512. But it is not! Does that change?' endif flatfile0=sxpar(guidhdr,'FLATFILE') if n_elements(flatfile) eq 0 then begin flatfile=flatfile0 endif else begin if flatfile0 ne flatfile then message,'Error: Guider images to be combined are using more than one flat frame. Are you sure the guider frames are all taken during one exposure? ' endelse darkfile0 = sxpar(guidhdr,'DARKFILE') if n_elements(darkfile_inmem) eq 0 then begin p0=strpos(darkfile0,'gimg') darkfile = strmid(darkfile0,p0,strlen(darkfile0)-p0) dark=mrdfits(guidedir+darkfile,0,/unsigned,hdr_dark) darkfile_inmem=darkfile0 endif else if darkfile0 ne darkfile_inmem then begin p0 = strpos(darkfile0,'gimg') darkfile = strmid(darkfile0,p0,strlen(darkfile0)-p0) dark=mrdfits(guidedir+darkfile,0,/unsigned,hdr_dark) darkfile_inmem=darkfile0 endif ;The guide camera has very little dark current. This is shown by comparing a 4 sec image with a 20 sec dark exposure. They are almost identical with a constant offset (the longer exposure image is lower in counts). But the bias has a shape: it gets stronger at small x. Therefore, we should use an offset dark image to bias-subtract the science image. ; Bias subtract the guider image darksize=size(dark,/dimen) darkx = darksize[0] & darky=darksize[1] if darkx eq imgsz[0] and darky eq imgsz[1] then begin biasdark=median(dark[darky:darkx-6*darky/512,*]) biasimg =median(img[darky:darkx-6*darky/512,*]) darksmooth=median(dark,dimen=2)#(fltarr(darky)+1) bias_to_sub = float(darksmooth)+(float(biasimg)-float(biasdark)) bimg = float(img)-bias_to_sub bimg = bimg-median(bimg) endif else begin splog,'Warning: Dark image size does not match the guider image size --- cannot use the dark image for dark subtraction.' splog,'Instead we will use simple median to subtract the image, which is less accurate, especially for fiber 7 and 3' bimg = float(img)-median(img) endelse if firstframe then begin totimg = bimg totexp = sxpar(guidhdr,'EXPTIME') endif else begin totimg=totimg+bimg totexp = totexp + sxpar(guidhdr,'EXPTIME') endelse firstframe = 0 if c2 eq 1 then goodproc=proc ; keep the name of the last proc-gimg file that exists. n_validframe +=1 endfor if n_validframe eq 0 then begin errorcode = 1 return endif info = mrdfits(guidedir+goodproc,6,/silent) if n_elements(flatfile_inmem) eq 0 then begin p0=strpos(flatfile0,'gimg') flatfile = strmid(flatfile0,p0,strlen(flatfile0)-p0) flat=mrdfits(guidedir+flatfile,0,/unsigned,hdr_flat) flatfile_inmem=flatfile0 endif else if flatfile0 ne flatfile_inmem then begin p0=strpos(flatfile0,'gimg') flatfile = strmid(flatfile0,p0,strlen(flatfile0)-p0) flat=mrdfits(guidedir+flatfile,0,/unsigned,hdr_flat) flatfile_inmem=flatfile0 endif ; Bias subtract the flat using the overscan if (size(flat,/dimen))[0] ne 1048 then begin message,'Error: I am expecting the size for the flat frame to be 1048 x1024. But it is not!' endif bflat = flat-median(flat[1024:1038,*]) ; Bin flat by 2x2 flatsize=size(bflat,/dimen) binflat = rebin(bflat,flatsize[0]/2,flatsize[1]/2.) thresh = ml_percentiles(binflat,96.5) bitmask = byte(binflat*0) bitmask[where(binflat le thresh)] = 4B ; binflat = binflat/median(binflat[where(bitmask eq 0B)]) binflat = binflat/ml_percentiles(binflat,99.8) ; flat-field the guider image fimg = totimg/binflat*(bitmask eq 0B)+totimg*(bitmask eq 4B) k=where(finite(fimg) eq 0,ct_inf) if ct_inf gt 0 then fimg[k]=0.0 ; readnoise 10.7 ADU, 17.1 e- errimg = sqrt(abs(fimg)+nimg*17.^2)*(bitmask eq 0B) + 1.e30*(bitmask eq 4B) if ct_inf gt 0 then errimg[k] = 1.e30 ; Make header for co-added guider image mkhdr,header,fimg sxdelpar,header,'DATE' sxaddpar,header,'TIMESYS',sxpar(sci_hdr,'TIMESYS'),before='COMMENT' sxaddpar,header,'TAI-BEG',sxpar(sci_hdr,'TAI-BEG'),before='COMMENT' sxaddpar,header,'DATE-OBS',sxpar(sci_hdr,'DATE-OBS'),before='COMMENT' sxaddpar,header,'GUIDER1',guider1,'The first guider image',before='COMMENT' sxaddpar,header,'GUIDERN',guidern,'The last guider image',before='COMMENT' sxaddpar,header,'TOTEXPTIME',totexp,'Total guider exposure time',before='COMMENT' sxaddpar,header,'OBSTIME',sxpar(sci_hdr,'EXPTIME'),'Total science exposure time',before='COMMENT' sxaddpar,header,'DARKFILE',darkfile0,'guider dark exposure',before='COMMENT' sxaddpar,header,'FLATFILE',flatfile0,'guider flat exposure',before='COMMENT' sxaddpar,header,'CARTID',sxpar(sci_hdr,'CARTID'),'The currently loaded cartridge',before='COMMENT' sxaddpar,header,'PLATEID',sxpar(sci_hdr,'PLATEID'),'The currently loaded plate',before='COMMENT' sxaddpar,header,'EXPOSURE',sxpar(sci_hdr,'EXPOSURE'),before='COMMENT' sxaddpar,header,'MJD',sxpar(sci_hdr,'MJD'),'APO fMJD day at start of exposure',before='COMMENT' sxaddpar,header,'DCNRA',dra,'applied user supplied offset in RA, arcsec',before='COMMENT' sxaddpar,header,'DCNDEC',ddec,'applied user supplied offset in Dec, arcsec',before='COMMENT' info = ml_trimtag(info,['XSTAR','YSTAR','DX','DY','DRA','DDec','FWHM','FLUX','MAG','SKY','SKYMAG','POSERR','STAMPSIZE','STAMPIDX']) ; Write the coadded guider file if ~keyword_set(dos) then begin ; check for directory, if doesn't exist, create it if file_test(outdir, /directory) eq 0 then spawn, '\mkdir -p '+outdir mwrfits,fimg,outdir+'cogimg-'+filestr+'.fits',header,/create mwrfits,errimg,outdir+'cogimg-'+filestr+'.fits' mwrfits,bitmask,outdir+'cogimg-'+filestr+'.fits' mwrfits,info,outdir+'cogimg-'+filestr+'.fits' endif else begin ; Do NOT write the coadded guider image ; For DOS - to update the output file header with information on seeing and transparency psfpar = psfparams(plateid, mjd, exposure, img=fimg, errimg=errimg, bitmask=bitmask, info=info, header=header, dos=dos, transparency=transparency) ; Use filename if provided if n_elements(filename) ne 0 then begin ; Check if filename exists res = file_test(filename) if res ne 1 then splog, 'ERROR: No DOS file found. Cannot add header information!' ; Grab header of output FITS file hdr = headfits(filename) sxaddpar, hdr, 'SEEING', strtrim(psfpar.fwhm,2), 'Seeing from Co-added Guider Image' sxaddpar, hdr, 'TRANSPCY', transparency, 'Transparency computed from Co-added Guider Image' modfits, filename, 0, hdr splog, 'Wrote seeing and transparency to output mgscisky DOS filename' endif else if n_elements(hdr) ne 0 then begin ; Have only the header info sxaddpar, hdr, 'SEEING', strtrim(psfpar.fwhm,2), 'Seeing from Co-added Guider Image' sxaddpar, hdr, 'TRANSPCY', transparency, 'Transparency computed from Co-added Guider Image' splog, 'Wrote seeing and transparency to output mgscisky DOS file' endif else splog, 'ERROR: No DOS filename or hdr provided. Skipping seeing and transparency output.' endelse return end