function read_vault,header,folder=folder,xrange=xrange,yrange=yrange,all=all,maps=maps,flatfile=flatfile ;+ ; NAME: ; READ_VAULT ; PURPOSE: ; Read VAULT2 data from a selected FoV ; EXPLANATION: ; Reads and extracts VAULT2 data from original level1 ; files into a variable, or map. ; OUTPUT data has 1+21 frames. First frame has a median value of all running frames ; ; CALLING SEQUENCE: ; data= read_vault([header,folder=folder,xrange=xrange,yrange=yrange,flatfile='flat.fits',/all,/map]) ; INPUTS: ; ; ; OPTIONAL INPUTS: ; FOLDER - Folder containing original Vault level1 files. ; [current folder is used if none is specified] ; XRANGE - Sun X range of subfield (in W-E Sun arcsec coordinates). ; Use only in combination with YRANGE. ; [Makes routine non-interactive] ; YRANGE - Sun Y range of subfield (in N-S Sun arcsec coordinates). ; Use only in combination with XRANGE. ; [Makes routine non-interactive] ; ; OPTIONAL OUTPUTS: ; HEADER - Variable to retreive header structure info ; ; OPTIONAL INPUT KEYWORD: ; ALL - Set this keyword to retreive all data in a single frame. ; To reduce output variable size, data is rescaled to UNIT data type (0 to 65535) ; ; MAPS - Set this keyword to return a D. ZARRO map as output. That is, an array of ; structures, each of which is a single map. ; ; FLAT - Correct flat fielding and vignetting, with a pre-computed estimation. ; EXAMPLE: ; Interactively select a subregion, and store as a variable and header info. ; (better for direct data manipulation) ; ; IDL> data=read_vault(header) ; ; Retreive all data in a single map variable. And plot map of meadian frame (first frame) ; (better for D. Zarro map routines) ; ; IDL> map=read_vault(/all,/maps) ; IDL> plot_map,map(0) ; ; Read an specified region (non interactive) and apply flat correction ; ; IDL> a=read_vault(flat='flat.fits',xrange=[-800,-60],yrange=[150,350]) ; ; MODIFICATION HISTORY: ; Written, Dec. 2008 Bruno S.-A.N. SSD / NRL ;- on_error,2 if ~keyword_set (folder) then folder='./' if ~keyword_set (maps) then if N_PARAMS() EQ 0 then message,'use /maps keyword or add variable to retrieve the header info, otherwise is not retreived: data=read_vault(header)',/cont if keyword_set(flatfile) then if file_test(flatfile) then begin flat=readfits(flatfile,/silent) endif else message,'flatfielding reqested but no flat file found' files=['composite.fits','Vault2_020614_181201UT.fits','Vault2_020614_181218UT.fits','Vault2_020614_181239UT.fits','Vault2_020614_181257UT.fits','Vault2_020614_181314UT.fits','Vault2_020614_181331UT.fits','Vault2_020614_181348UT.fits','Vault2_020614_181405UT.fits','Vault2_020614_181422UT.fits','Vault2_020614_181439UT.fits','Vault2_020614_181456UT.fits','Vault2_020614_181513UT.fits','Vault2_020614_181530UT.fits','Vault2_020614_181547UT.fits','Vault2_020614_181604UT.fits','Vault2_020614_181621UT.fits','Vault2_020614_181639UT.fits','Vault2_020614_181656UT.fits','Vault2_020614_181713UT.fits','Vault2_020614_181730UT.fits','Vault2_020614_181747UT.fits'] message,'checking if all '+strtrim(string(n_elements(files)),2)+' needed files are there',/cont for i=0,n_elements(files)-1 do if ~file_test(folder+files[i]) then message,'File missing: '+files[i] message,'Reading composite frame',/cont T = SYSTIME(1) fits2map,folder+'composite.fits',map DT=SYSTIME(1) - T ;restore,folder+'/poster.save' ;map=c map.data=sigrange(map.data) ;select subframe if keyword_set (xrange) AND keyword_set (yrange) then sub_map,map,smap,xrange=xrange,yrange=yrange $ else if keyword_set (all) then smap=map $ else begin device,decomposed=0 window,0,xs=460*3,ys=300*3 sub_map,map,smap,grid=7,/plot,title='Select region' endelse si=size(smap.data) ;announce size of variable mem=n_elements(smap.data)/1048576. if keyword_set (all) then mem*=16 else mem*=32 message,'output will use up to'+strcompress(fix(mem))+' Mb of memory',/cont si=size(smap.data) ;initilize data output if keyword_set (all) then begin data=uintarr(si[1],si[2],22) data[*,*,0]=uint(smap.data >0 ) endif if keyword_set(maps) then data=smap else $ begin data=fltarr(si[1],si[2],22) data[*,*,0]=smap.data endelse vault={xc:smap.xc,yc:smap.yc,dx:smap.dx,dy:smap.dy,time:replicate(smap.time,22),id:smap.id,dur:replicate(smap.dur,22),xunits:smap.xunits,yunits:smap.yunits,roll_angle:smap.roll_angle,roll_center:smap.roll_center,l0:smap.l0,b0:smap.b0,Rsun:smap.Rsun} vault.time[0]='median frame' vault.dur[0]=1 print,'% READING the selected 21 subfields.' DT2=0 for i=1,21 do begin writeu,-1,strtrim(string(i),2)+',' T2 = SYSTIME(1) mreadfits,folder+files(i),a,b,/silent DT2+=SYSTIME(1) - T if keyword_set(flatfile) then begin message,'flatfielding data...',/cont bs=size(b)-1 fs=size(flat)-1 ;flat frame can have few more or few less lines or col than frame flatt=fltarr(bs[1]+1,bs[2]+1) flatt[0:bs[1]0,bs[2]-fs[2]>0) ;find shifts between data and flat, if any sh=shc(flatt[1300:1500,1300:1500],b[1300:1500,1300:1500] LT 0.3*mean(b[1300:1500,1300:1500])) ;print,sh flatt=shift(flatt,-sh) bf=b/flatt flatt=0 ;remove points NaN or divided by 0 b=bfmin(b)/2 ;remove NaNs and 1/0s endif index2map,a,b,imap sub_map,imap,smapi,ref_map=smap,/noplot,err=err IF err EQ 'Insufficient points in sub-region' then message,'skip frame '+strtrim(string(i),2),/cont else begin imap=coreg_map(smapi,smap) tmp=imap.data if keyword_set (all) then data[*,*,i]=uint(tmp/max(tmp)*65535 >0 ) if keyword_set(maps) then data=merge_struct(data,imap) else data[*,*,i]=imap.data vault.time[i]=imap.time vault.dur[i]=imap.dur plot_map,smapi,grid=7 endelse endfor print,'.' ;print,'Actual time: '+strtrim(string(fix(DT2-DT)),2)+' s.' header=vault if keyword_set (all) then message,'Warning: Data output is a rescaled UINT array',/cont ;convert data to array of map structures return,data end