; file: sdo_stx_combine_v1.pro ; init: Jan 23 2018 Rob Rutten Deil ; last: Apr 4 2018 Rob Rutten Deil ; note: v1 was original for combining DST/IBIS and DST/MXIS ;+ pro sdo_stx_combine,indirs,selnames,outdir,$ outcadence=outcadence,outduration=outduration,$ outpixel=outpixel,outnx=outnx,outny=outny,$ splinip=splinip ; PURPOSE: ; combine STX2SDO datasets, each SDO-aligned, into one co_set ; i.e., synchronized and co-spatial = ready for showex.pro ; also add corresponding SDO files ; DESCRIPTION: ; distill input timing, pixel size, dimensions from header ; specify output timing, pixel size, dimensions ; rewrite all STX cubefiles together into outdir ; ?? cadence = A, B, fastest, mean? ; ?? duration = A, B, overlap? ; ?? pizel size = A, B, smallest, largest? ; ?? dimension extent = A, B, smallest, largest? ; ; CALL: ; see above ; ; INPUTS: ; indirs: array of strings with directory paths ; selnames: array of strings with fits file-selector substrings ; outdir: string with directory path ; ; OPTIONAL KEYWORD INPUTS: ; not yet there; now smallest px, fastest cadence, common area ; ; OUTPUTS: ; files in outdir ; ; HISTORY: ; Jan 23 2018 RR: start = IBIS and MXIS on SDO ; Mar 10 2018 RR: more than 2 STX ;- ; answer no-parameter query if (n_params(0) lt 3) then begin sp,sdo_stx_combine return endif ; set keyword defaults if (n_elements(outcadence) eq 0) then outcadence=-1 if (n_elements(outduration) eq 0) then outduration=-1 if (n_elements(outpixel) eq 0) then outpixel=-1 if (n_elements(outnx) eq 0) then outnx=-1 if (n_elements(outny) eq 0) then outny=-1 if (n_elements(splinip) eq 0) then splinip=0 ; set wall-clock timer (seconds) timestart=systime(1) ; print pacifier print,' ===== sdo_stx_combine started at UT = '+$ anytim2utc(timestart,/vms) ; numbers ndirs=n_elements(indirs) nselnames=n_elements(selnames) ; make outdir if not yet existing spawn,'mkdir -p '+outdir ; define keyword arrays bitpix_arr=intarr(ndirs) nx_arr=intarr(ndirs) ny_arr=intarr(ndirs) nt_arr=intarr(ndirs) starttai_arr=dblarr(ndirs) endtai_arr=dblarr(ndirs) cadence_arr=fltarr(ndirs) cdelt1_arr=fltarr(ndirs) cdelt2_arr=fltarr(ndirs) xcen_arr=fltarr(ndirs) ycen_arr=fltarr(ndirs) px_arr=fltarr(ndirs) ; loop over indirs to get keyword parameters from first fits file in there for idir=0,ndirs-1 do begin files=findfile(indirs[idir]+'/*.fits') nfiles=n_elements(files) if (nfiles eq 0) then begin print,' ##### sdo_stx_combine abort: no files in indir = '+indirs[idir] return endif header=headfits_rr(files[0]) ; get this dir parameters from first-file header in this dir header=headfits_rr(files[0]) bitpix_arr[idir]=fxpar(header,'bitpix') nx_arr[idir]=fxpar(header,'naxis1') ny_arr[idir]=fxpar(header,'naxis2') nt_arr[idir]=fxpar(header,'naxis3') starttai_arr[idir]=anytim2tai(fxpar(header,'starttim')) cadence_arr[idir]=fxpar(header,'cadence') cdelt1_arr[idir]=fxpar(header,'cdelt1') cdelt2_arr[idir]=fxpar(header,'cdelt2') xcen_arr[idir]=fxpar(header,'xcen') ycen_arr[idir]=fxpar(header,'ycen') ; pixels px_arr[idir]=(cdelt1_arr[idir]+cdelt2_arr[idir])/2. ; end times endtai_arr[idir]=starttai_arr[idir]+(nt_arr[idir]-1)*cadence_arr[idir] endfor ; end loop over indirs to get parameters ; define spatial output parameters ; default = overlap around shifted-together centers, finest pixel scale px_out=min(px_arr) pxratio_arr=px_out/px_arr nxcut_arr=nx_arr/pxratio_arr nycut_arr=ny_arr/pxratio_arr nxcut_out=fix(min(nxcut_arr)) nycut_out=fix(min(nycut_arr)) ; define output timing ; default = all-facility overlap period at fastest cadence starttai_out=max(starttai_arr) endtai_out=min(endtai_arr) cadence_out=min(cadence_arr) nt_out=fix((endtai_out-starttai_out)/cadence_out) timearr_out=starttai_out+indgen(nt_out)*cadence_out ; measure shifts for one pair. Shift must be measured because stx2sdo ; recenters SDO to have each STX field centered in the stx2sdo field. ; Time gets the same by using timearr_out. ; find SDO hmicont files for shift determination shiftfiles=strarr(ndirs) for idir=0,ndirs-1 do begin files=findfile(indirs[idir]+'/*.fits') nfiles=n_elements(files) for ifile=0,nfiles-1 do $ if (strmatch(files[ifile],'*hmicont*')) then shiftfiles[idir]=files[ifile] if (shiftfiles[idir] eq '') then begin print,' ##### sdo_stx_combine abort: no hmicont file in dir'+indirs[idir] return endif endfor ; cut shiftfiles to short duration to find shifts shortfiles='/tmp/shortfile_'+trim(string(indgen(ndirs)))+'.fits' timearr_short=timearr_out[0:4] for idir=0,ndirs-1 do begin timearr_in=starttai_arr[idir]+indgen(nt_arr[idir])*cadence_arr[idir] reformcubefile,shiftfiles[idir],shortfiles[idir],verbose=0,$ congridfactor=1./[pxratio_arr[idir],pxratio_arr[idir]],$ cutcentralx=nxcut_out,cutcentraly=nycut_out,$ timearr_in=timearr_in,timearr_out=timearr_short endfor ;; ; inspect ;; showex,shortfiles ; indeed sizable shifts ;; STOP ; find the shifts from first images, pairwise dirs to first dir file dirshifts=fltarr(ndirs,2) for idir=1,ndirs-1 do begin im_A=readfitscube(shortfiles[0],trange=[0,0]) im_B=readfitscube(shortfiles[idir],trange=[0,0]) shift_AB=findimshift_rr(im_A,im_B,/subpix,filter=10,cropborders=1) print,' ===== '+shortfiles[idir]+' has shift to '+shortfiles[0]+' = '+$ trimd(shift_AB) dirshifts[idir,0]=shift_AB[0] dirshifts[idir,1]=shift_AB[1] endfor ; ----- process all files for idir=0,ndirs-1 do begin timearr_in=starttai_arr[idir]+indgen(nt_arr[idir])*cadence_arr[idir] ; in first dir do all .fits files including all SDO files if (idir eq 0) then begin inselstring='.fits' outtail='.fits' doallfiles,indirs[idir],outdir,inselstring,outtail,$ 'reformcubefile',$ congridfactor=1./[pxratio_arr[idir],pxratio_arr[idir]],$ shift=[dirshifts[idir,0],dirshifts[idir,1]],$ cutcentralx=nxcut_out,cutcentraly=nycut_out,$ timearr_in=timearr_in,timearr_out=timearr_out,$ missingvalue=-1,splinip=splinip endif else begin ; in other dirs do only the selnames .fits files (to exclude SDO) for isel=0,nselnames-1 do begin inselstring='*'+selnames[isel]+'*.fits' outtail='.fits' doallfiles,indirs[idir],outdir,inselstring,outtail,$ 'reformcubefile',$ congridfactor=1./[pxratio_arr[idir],pxratio_arr[idir]],$ shift=[dirshifts[idir,0],dirshifts[idir,1]],$ cutcentralx=nxcut_out,cutcentraly=nycut_out,$ timearr_in=timearr_in,timearr_out=timearr_out,$ missingvalue=-1,splinip=splinip endfor endelse ; do all La Palma .icube files in every dir ; (keep them .icube to have them recognized by movex_loadfiles.pro) inselstring='.icube' outtail='.icube' doallfiles,indirs[idir],outdir,inselstring,outtail,$ 'reformcubefile',$ congridfactor=1./[pxratio_arr[idir],pxratio_arr[idir]],$ shift=[dirshifts[idir,0],dirshifts[idir,1]],$ cutcentralx=nxcut_out,cutcentraly=nycut_out,$ timearr_in=timearr_in,timearr_out=timearr_out,$ missingvalue=-1,splinip=splinip endfor ; end loop over all dirs ; print elapsed time timedone=systime(1) timelaps=ntostr((timedone-timestart)/60.,format='(F11.1)') print,' ===== sdo_stx_combine took '+timelaps+' minutes' ; print cadence (for time delay setting in showex) print,' ----- cadence_out ='+trimd(cadence_out) end ; =============== main for testing per IDLWAVE H-c ====================== ; Schad DST 10830 data ;; cd,'/media/rutten/RRDATA/alldata/DST/2017-06-14-schad' ;; indirs=['ibisB2sdo','mxis2sdo','iris2sdo_1400'] ;; selnames=['ibis','mxis','iris'] ;; outdir='ibis_mxis_iris_sdo/' ;; splinip=1 ; Luc SST contrail data now put on rotated SDO cd,'/media/rutten/SSTDATA/alldata/SST/2014-06-21-quiet' indirs=['sst2sdo','iris2sdo'] selnames=['crispex','wb','iris'] outdir='sst_iris_sdo' splinip=1 sdo_stx_combine,indirs,selnames,outdir,splinip=splinip end