pro eit_prep,filename,outheader,outimage, $ ; ;----------------------Calibration, processing options no_calibrate = no_calibrate, cosmic=cosmic,$ surround=surround,fill_block=fill_block,$ no_index=no_index,$ ; ;----------------------Output controls outfits=outfits,nodata=nodata,outdir=outdir,$ ; ;----------------------Misc. verbose=verbose,save_zero=save_zero, $ nrl=nrl,n_block=n_block,data=data, $ image_no=image_no, noprompt = noprompt,no_roll=no_roll, _extra=extra ;+ ; NAME: ; EIT_PREP ; ; PURPOSE: ; Process an EIT image. ; ; Processing routine for preparing Level-1 EIT data. ; Presently this version reads in the named FITS ; image and produces a background subtracted, degridded, ; flat-fielded, and degradation corrected output array or FITS file. ; In the future this routine will also vignette correct. ; The input may also be the output index structure and data cube as ; produced by read_eit. ; ; CATEGORY: ; FITS processing ; ; CALLING SEQUENCE: ; EIT_PREP,filename,[outheader,outimage] [data=data] [,/outfits] ; [,/verbose] [,/cosmic] [,/no_calibrate] [,/no_index] ; [,/save_zero] [,/nrl] [n_block = n_block] ; [,/surround] [,/nodata] [,outdir = outdir] ; [,image_no = image_no] ; ; INPUTS: ; filename - FITS file name or the output of eit_catrd ; may also be index structure (e.g. output of read_eit) ; ; ; OPTIONAL INPUT KEYWORD PARAMETERS: ; outfits - set if want to produce an output FITS file ; (or set string filename) ; outdir - Destination directory for prepped FITS files ; (if outdir is set, /outfits is assumed) ; nrl - set if input file was processed at NRL ; surround - set to replace missing block data with the average of ; surrounding blocks. ; image_no - Image number(s) (or 'ALL') to process in multiple image LZ file. ; fill_block - intarr(32,32) which is used to replace missing data blocks ; save_zero- set to retain images that have only 0 counts. ; verbose - set for a few more messages ; data - data array(cube) corresponding to index structure ; nodata - set if don't want output array (auto set if lt 3 params) ; cosmic - set if want to remove cosmic rays - may remove small ; real features ; no_calibrate - set if want "raw" returned images, i.e. only background subtracted ; no_index - set if don't want to use recommended Index correlation ; for final response scaling ; no_roll - do not correct for SOHO roll (mainly for testing) ; ; OPTIONAL OUTPUTS: ; outheader- FITS header or output index structure ; outimage - Processed EIT image. Default data type is Float ; ; OPTIONAL OUTPUT KEYWORD PARAMETERS: ; n_block - Number of blocks repaired ; ; COMMON BLOCKS: none. ; ; SIDE EFFECTS: ; If /outfits is set, then a FITS file is written. ; ; RESTRICTIONS: ; 1) Processes only one 3-Dim file at a time, either a single image ; vector, or the entire file. ; ; 2) Does not process a mixture of 2D and 3D files ; ; 3) If filename is a vector, there is no checks to make sure ; that all the images are the same size. User must be careful ; to check this. ; ; PROCEDURE: ; The FITS file is read in. Any missing blocks are replaced ; with the present dark noise value. The dark noise is ; subtracted off the entire image. The image is degridded. ; The image is flat fielded and corrected for degradation. ; The image is optionally returned. ; A new FITS file with an updated header is written out ; with the format eit_l1_yyyymmdd_hhmmss.fits ; ; Note: the De-gridding and Flat-fielding are multiplicative ; operations with an AVERAGE value over the full field ; of 1.0 ; ; MODIFICATION HISTORY: ; Written by: J. Newmark Date. July 1996 ; 1996 Sep 24 - J. Newmark - add flat fielding ; 1996 Oct 8 - J. Newmark - allow NRL format files ; 1996 Oct 16, J. R. Lemen, Changed calling arguments. ; 1996 Dec 09 - J. Newmark - add usage of image_no throughout. ; 1997 Feb 10 - J. Newmark - allow input if index structure and data ; cube ; 1997 Mar 09 - J. Newmark - write out FITS files from indx/data ; structures ; 1997 Mar 10 - J. Newmark - fix output FITS file names. ; 1997 Apr 09 - J. Newmark - account for new exptime ; 1997 Apr 12 - J. Newmark - account for correction of Al+1 to ClEAR ; 1997 Aug 14 - J. Newmark - implemented Handy requests of outdir, use ; of nodata ; 1997 Aug 18 - J. Newmark - updated conversion of struct->fitshead ; fix up code for /surround ; clean up code for output file name ; 1997 Oct 06 - J. Newmark - added preliminary response normalization ; 1997 Oct 30 - J. Newmark - update header for LZ Partial Fov images ; 1997 Nov 10 - J. Newmark - add new grids after July 1997 ; ; Version 2.0 ; 1998 Mar 27 - J. Newmark - support ALL or vector images in single 3D ; file, speed up 3D processing, ; clean up history, pathway checks, ; update documentation, fix n_block bug ; added adjust keyword ; 1998 Jun 3 - D. Zarro (SAC/GSFC) ; - sprinkled liberally with calls to temporary ; 1998 Jul 22 - J. Newmark - add phot_norm keyword ; 1998 Jul 23 - J. Newmark - add cosmic keyword ; 1998 Dec 08 - J. Newmark - change adjust keyword ; 2000 Jan 05 - J. Newmark - new response normalization, deleted ; phot_norm keyword ; 2000 Aug 04 - J. Newmark - NOPROMPT keyword for batch processing ; large data sets, cleaned up memory managment ; 2001 Mar 15 - F. Auchere - added time varying dark current ; ;Version 3.0 ; 2001 May 18 - J. Newmark/F. Auchere - new response normalization ; 2001 Jun 21 - N. Rich - Add NRL check at (near) line 332 ; 2001 Aug 15 - Zarro (EITI/GSFC) - add optional FITS out file name in ; OUTFITS ; ;Version 4.0 ; 2001 Dec 10 - J. Newmark - changed default options, ; new calibration, new default options lead ; to default output type as float ; 2002 Jan 12 - F. Auchere - modified call to eit_dark to handle ; different offset in binned images ; 2003 Apr 17 - D. Biesecker - modified call to is_fits to correctly ; handle collapsed directory path names ; 2003 Jul 17 - C.A. Young - added 180 degree roll correction ; 2003 Sep 13 - C.A. Young - added fix for multifile list ; 2004 Feb 4 - D.M. Zarro - Update header for rotation correction ;- utc_date_19960327 = anytim2utc('1996/03/27') utc_date_19960323 = anytim2utc('1996/03/23') utc_obe_good = anytim2utc('1996/07/18') t_obe = utc_obe_good.mjd + 1.d-3*utc_obe_good.time/86400. flag_3d = 0 use_fits = 1 filter_factor = [0.49, 0.49, 0.33, 0.29] wavelength = [171, 195, 284, 304] if is_string(outfits) then outfile=outfits if n_elements(outdir) then outfits = 1 else outdir = '' yes_data = not keyword_set(nodata) or n_params() eq 3 t0 = systime(1) ; Keep track of running time if datatype(filename) eq 'STC' then begin if keyword_set(data) then begin use_fits = 0 image_no = 0 endif else begin message, /info, 'No Data supplied, reading FITS files...' filename = filename.filename endelse endif nfiles = n_elements(filename) vv = fltarr(nfiles) n_block = intarr(nfiles) if use_fits then begin stat = is_fits(filename(0)) if stat eq 0 then begin files = eit_file2path(filename) stat = is_fits(files(0)) if stat eq 0 then begin files = eit_file2path(filename,/collapse) stat= is_fits(files(0)) if stat eq 0 then begin message,'The files do not appear to be FITS file.',/info return endif endif endif else files = filename if yes_data then begin noprompt = keyword_set(noprompt) or (get_logenv('ssw_batch') ne '') prompt = 1 - noprompt tsize = total(file_stat(files,/size))/1.e6 if tsize ge 100 and prompt then begin print,'You have requested '+string(tsize)+' MB of data' print,'You can select the output to be only FITS files if this '+$ 'is too large for your machine"s memory.' ans='' read,ans,prompt='Enter C (continue) or F (FITS only): ' if strupcase(ans) eq 'F' then begin yes_data = 0 outfits = 1 endif endif endif if keyword_set(image_no) and nfiles gt 1 then begin message,'This program only works on 1 3-Dim file at a time.',/info return endif endif i = 0 repeat begin ; MAIN LOOP: begins if ((i+1) mod 10) eq 0 then message,/info, $ string('Working on',i+1,'/',strtrim(nfiles,2), $ ' ETC=',(nfiles-i+1)*(systime(1)-t0)/60./float(i+1), $ ' mins.',form='(a,i4,3a,f5.1,a)') if use_fits then begin thisfile = files(i) if n_elements(image_no) eq 0 then begin image = readfits(thisfile,header,silent=1-keyword_set(verbose)) image_no = 0 endif else begin if flag_3d then image_no = images3d(i) header = headfits(thisfile) if n_elements(image_no) gt 1 then begin nfiles = n_elements(image_no) files = replicate(files,nfiles) images3d = image_no image_no = images3d(0) flag_3d = 1 vv = fltarr(nfiles) n_block = intarr(nfiles) endif else if strupcase(image_no) eq 'ALL' then begin nfiles = eit_fxpar(header,'naxis3') files = replicate(files,nfiles) images3d = indgen(nfiles) image_no = 0 flag_3d = 1 vv = fltarr(nfiles) n_block = intarr(nfiles) endif if image_no eq 0 then begin image = readfits(thisfile,silent=1-keyword_set(verbose)) endif else begin image = readfits(thisfile,silent=1-keyword_set(verbose),$ nslice = image_no) endelse endelse if keyword_set(nrl) then header = nrl2eit(header) endif else begin header = filename(i) thisfile = eit_file2path(header.filename) image = data(*,*,i) endelse ; set final keyword for degrid routine if strpos(strlowcase(thisfile), 'efz') ge 0 or strtrim(eit_fxpar(header,$ 'DATASRC'),2) eq 'LZ file' then final = 1 else final = 0 date_obs = eit_fxpar(header,'date_obs',image_no=image_no) ;== Added 15 March 2001 F. Auchere utc_date_obs = anytim2utc(date_obs) hist = eit_fxpar(header,'history') IF keyword_set(NRL) THEN version = 3.0 ELSE version = float(strmid(hist(0),8,3)) ;==Add NRL check, 21Jun01, nbr ; Fix missing blocks, dark subtract: if keyword_set(fill_block) and keyword_set(surround) then $ message,'Cannot set both fill_bock and surround' if keyword_set(cosmic) then begin image = cosmicr(image) tagval = 'Replaced Cosmic Rays with Median Filtering' if use_fits then fxaddpar, header,'HISTORY', tagval else $ header = rep_tag_value(header,[header.history,tagval],'HISTORY') endif bin = round(eit_fxpar(header, 'CDELT1')/eit_pixsize()) case 1 of keyword_set(surround): begin if not keyword_set(no_degrid) then $ image = raise_missing_blocks(image, repl=eit_dark(date_obs,bin=bin),n_block=nblock,/no_copy) $ else image = raise_missing_blocks(image, /surround,n_block=nblock,/no_copy) ;== Modified 15 March 2001 F. Auchere tagval = 'Replaced missing blocks with SURROUND' end keyword_set(fill_block) : begin image =raise_missing_blocks(image,fill_block=fill_block,n_block=nblock,/no_copy) tagval = 'Replaced missing blocks with FILL_BLOCK' end else: begin image = raise_missing_blocks(image, repl=eit_dark(date_obs, bin=bin),n_block=nblock,/no_copy) ;== Modified 15 March 2001 12 jan 2002 F. Auchere tagval = 'Replaced missing blocks with EIT_DARK' end endcase n_block(i) = nblock if use_fits then fxaddpar, header,'HISTORY', tagval else $ header = rep_tag_value(header,[header.history,tagval],'HISTORY') image = temporary(image) - eit_dark(date_obs, bin=bin) ;== Modified 15 March 2001, 12-jan 2002 F. Auchere tagval = 'Subtracted dark current, EIT_DARK' if use_fits then fxaddpar, header,'HISTORY', tagval else $ header = rep_tag_value(header,[header.history,tagval],'HISTORY') ; calibrate - degrid, flat field, normalize by exposure time, time varying ; response correction if not keyword_set(no_calibrate) then begin image = eit_degrid(image, header, final = final, image_no = image_no,$ /no_copy) if keyword_set(surround) then $ image = raise_missing_blocks(image, /surround,n_block=nblock,/no_copy) tagval = 'Degridded image with EIT_DEGRID' if use_fits then fxaddpar, header,'HISTORY', tagval else $ header = rep_tag_value(header,[header.history,tagval],'HISTORY') image = eit_flat(image, header, image_no = image_no,/no_copy) tagval = 'Flat fielded image with EIT_FLAT' if use_fits then fxaddpar, header,'HISTORY', tagval else $ header = rep_tag_value(header,[header.history,tagval],'HISTORY') if version(0) ge 2.0 then shut_time = 0 else begin t_obs = utc_date_obs.mjd + 1.d-3*utc_date_obs.time/86400 shut_time = eit_fxpar(header,'SHUTTER CLOSE TIME',image_no=image_no) if shut_time le 0 or t_obs lt t_obe then shut_time = 2.1 endelse exptime=eit_fxpar(header,'exptime',image_no=image_no)+shut_time image = temporary(image) / exptime tagval = 'Exposure normalized (per sec)' if use_fits then fxaddpar, header,'HISTORY', tagval else $ header = rep_tag_value(header,[header.history,tagval],'HISTORY') filter = strupcase(strtrim(eit_fxpar(header,'filter',image_no=image_no),2)) wave = eit_fxpar(header,'wavelnth',image_no=image_no) i_wave = where(wavelength eq wave) if (not final) and (utc_date_obs.mjd lt utc_date_19960327.mjd) and $ (utc_date_obs.mjd ne utc_date_19960323.mjd) then begin if filter eq 'CLEAR' then image = temporary(image) / filter_factor(i_wave(0)) endif else if filter eq 'AL +1' then image = temporary(image)/filter_factor(i_wave(0)) if n_elements(leak_284) ne 0 and filter eq 'CLEAR' and wave eq 284 then $ image = temporary(image)/leak_284 tagval = 'Flux normalized to CLEAR Filter' if use_fits then fxaddpar, header,'HISTORY', tagval else $ header = rep_tag_value(header,[header.history,tagval],'HISTORY') wave = eit_fxpar(header,'wavelnth',image_no=image_no) image = temporary(image) / eit_norm_response(utc_date_obs,wave,header,$ no_index=no_index) tagval = 'Flux normalized for time response changes' if use_fits then fxaddpar, header,'HISTORY', tagval else $ header = rep_tag_value(header,[header.history,tagval],'HISTORY') endif else image = fix(round(temporary(image))) ; ; FITS header updates, needed for LZ partial field of view. if use_fits and final and eit_fxpar(header,'naxis3') gt 0 then begin fxaddpar, header, 'WAVELNTH', eit_fxpar(header,'wavelnth',image_no=image_no) fxaddpar, header, 'CFTEMP', eit_fxpar(header,'cftemp',image_no=image_no) fxaddpar, header, 'DATE_OBS', eit_fxpar(header,'date_obs',image_no=image_no) fxaddpar, header, 'FILTER', eit_fxpar(header,'filter',image_no=image_no) fxaddpar, header, 'EXPTIME', eit_fxpar(header,'exptime',image_no=image_no) up_corr = where(strpos(header,'COMMENT CORREC') ne -1) if up_corr(0) ne -1 then header(up_corr) = 'COMMENT '+$ 'CORRECTED DATE_OBS = '+ anytim2utc(eit_fxpar(header,$ 'corrected date_obs',image_no=image_no),/ccsds) endif ; Pointing adjustment - 3/30/98 JSN if version lt 2.6 and use_fits then begin fxaddpar, header, 'CDELT1', eit_pixsize() fxaddpar, header, 'CDELT2', eit_pixsize() xy = eit_point(utc_date_obs) fxaddpar, header, 'CRPIX1', xy(0) fxaddpar, header, 'CRPIX2', xy(1) endif ; Set up the output file name: dat_ob=anytim2cal(eit_fxpar(header,'date_obs',image_no=image_no),form=8) if is_string(outfile) then outname=outfile else $ outname = 'eit_l1_'+strmid(dat_ob,0,8)+'_'+strmid(dat_ob,8,6)+'.fits' outname = concat_dir(outdir,outname) ; Write the FITS file if keyword_set(outfits) then begin if use_fits then writefits, outname, image, header else begin hdr = struct2fitshead(header) writefits,outname,image,hdr delvarx,hdr endelse endif ; Concatenate the headers and images if i eq 0 then outheader = header else begin if use_fits then begin szo = size(outheader) szh = size(header) case 1 of szh(1) lt szo(1): header = [header,' '] szh(1) gt szo(1): header = header(0:szo(1)-1) else: endcase outheader = [[outheader],[header]] endif else begin out_str = outheader(0) outheader = [str_copy_tags(out_str,outheader), str_copy_tags(out_str,header)] endelse endelse if i eq 0 and yes_data then begin ss = size(image) if keyword_set(no_calibrate) then outimage = intarr(ss(1),ss(2),nfiles) else $ outimage = fltarr(ss(1),ss(2),nfiles) endif vv(i) = max(image) if yes_data then outimage(0,0,i) = temporary(image) i = i + 1 endrep until i eq nfiles ; MAIN LOOP ends if not keyword_set(save_zero) then begin ; Reject blank images jj = where(vv ne 0,ncount) if ncount eq 0 then begin message,/info,'All images are zero' outimage = 0 outheader = '' endif else begin if yes_data and ncount ne nfiles then outimage = outimage(*,*,jj) if use_fits then outheader = outheader(*,jj) else $ outheader = outheader(jj) n_block = n_block(jj) if n_elements(jj) eq 1 then n_block = n_block(0) endelse endif ;-- correct for 180 degree roll ;-- check in case OUTHEADER is a structure. If so, we handle it ; differently. if (1-keyword_set(no_roll)) then begin in_struct=have_tag(outheader,'SC_ROLL',rindex) for i=0,nfiles-1 do begin if in_struct then begin header=outheader[i] sc_roll=header.(rindex) endif else begin header=outheader[*,i] sc_roll = fxpar(header, 'SC_ROLL') endelse if abs(fix(sc_roll)) eq 180 then begin outimage[*,*,i]=rotate(outimage[*,*,i],2) outheader[*,i]=rot_fits_head(header) endif end endif message,/info,string('took',systime(1)-t0,' seconds',format='(a,f5.1,a)') end
Last revised: - Wed May 9 21:45:02 2007- F. Auchère