; file: sdo_stx_combine_v2.pro ; init: Jan 23 2018 Rob Rutten Deil ; last: Jan 19 2020 Rob Rutten Deil ; note: v2 = SST CRISP+CHROMIS combining separately stx2rotsdo aligned files ;+ pro sdo_stx_combine,stxnames,outdir,$ selectfiles=selectfiles,granfiles=granfiles, $ outcadence=outcadence,outduration=outduration,$ outpixel=outpixel,outnx=outnx,outny=outny,$ splinip=splinip ; PURPOSE: ; combine STX2ROTSDO data sequences, 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: ; stxnames: array STX codes defining stx2rotsdo dirs and sdo_stx_align dirs ; outdir: string with directory path ; ; OPTIONAL KEYWORD INPUTS: ; selectfiles = string array path/file names: do only these ; granfiles = STX showing granulation, default *wb*, *cont" ; others not yet done; should define above ?? output choices. ; now uses smallest px, fastest cadence, largest area ; ; OUTPUTS: ; files in outdir ; ; HISTORY: ; Jan 23 2018 RR: v1 start on IBIS + MXIS ; Mar 10 2018 RR: v1 more than 2 STX ; Jan 12 2020 RR: v2 = combine stx2rotsdo resuls ;- ; answer no-parameter query if (n_params(0) lt 2) then begin sp,sdo_stx_combine return endif ; set keyword defaults if (n_elements(selectfiles) eq 0) then selectfiles='' if (n_elements(granfiles) eq 0) then granfiles='' 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 nstx=n_elements(stxnames) ; make outdir if not yet existing spawn,'mkdir -p '+outdir ; define keyword arrays ;RR ?? no angle? bitpix_arr=intarr(nstx) nx_arr=intarr(nstx) ny_arr=intarr(nstx) nt_arr=intarr(nstx) starttai_arr=dblarr(nstx) endtai_arr=dblarr(nstx) cadence_arr=fltarr(nstx) cdelt1_arr=fltarr(nstx) cdelt2_arr=fltarr(nstx) xcen_arr=fltarr(nstx) ycen_arr=fltarr(nstx) px_arr=fltarr(nstx) angle_arr=fltarr(nstx) ; loop over instx to get keyword parameters from hmicont file in each for istx=0,nstx-1 do begin files=findfile(stxnames[istx]+'2rotsdo/hmicont*.fits') nfiles=n_elements(files) if (nfiles eq 0) then begin print,' ##### sdo_stx_combine abort: no hmicont file in indir = '$ +instx[istx] return endif header=headfits_rr(files[0]) ; get this dir parameters from hmicont header in this dir header=headfits_rr(files[0]) bitpix_arr[istx]=fxpar(header,'bitpix') nx_arr[istx]=fxpar(header,'naxis1') ny_arr[istx]=fxpar(header,'naxis2') nt_arr[istx]=fxpar(header,'naxis3') starttai_arr[istx]=anytim2tai(fxpar(header,'starttim')) cadence_arr[istx]=fxpar(header,'cadence') cdelt1_arr[istx]=fxpar(header,'cdelt1') cdelt2_arr[istx]=fxpar(header,'cdelt2') xcen_arr[istx]=fxpar(header,'xcen') ; unreliable if rotated ycen_arr[istx]=fxpar(header,'ycen') ; unreliable if rotated angle_arr[istx]=fxpar(header,'crota2') ; pixels px_arr[istx]=(cdelt1_arr[istx]+cdelt2_arr[istx])/2. ; end times endtai_arr[istx]=starttai_arr[istx]+(nt_arr[istx]-1)*cadence_arr[istx] endfor ; end loop over stxnames to get parameters ; define spatial output parameters ; default = overlap around shifted-together centers, finest pixel scale px_out=min(px_arr) ; finest pixel size pxratio_arr=px_out/px_arr nxcut_arr=nx_arr/pxratio_arr nycut_arr=ny_arr/pxratio_arr nxcut_out=fix(min(nxcut_arr)) ; smallest FOV size (largest > error) nycut_out=fix(min(nycut_arr)) ; define output angle = mean of stx2rotsdo values angle_out=avg(angle_arr) delangle_arr=angle_arr-angle_out ; define output timing ; default = all-sequence overlap period at fastest cadence starttai_out=max(starttai_arr) ; latest start endtai_out=min(endtai_arr) ; earliest end cadence_out=min(cadence_arr) ; fastest cadence 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. ; Times gets the same by using timearr_out. Tricky Dicky since needed ; precision is way below SDO resolution. In instantaneous comparisons ; the individual synchronized SDO images show considerable granular ; scale morphology variations likely from all the interpolations ; (already at JSOC for the cutouts). Smear and time-averaging proved ; necessary. ; find SDO hmicont files for shift determination shiftfiles=strarr(nstx) for istx=0,nstx-1 do begin files=findfile(stxnames[istx]+'2rotsdo/*.fits') nfiles=n_elements(files) for ifile=0,nfiles-1 do $ if (strmatch(files[ifile],'*hmicont*')) then shiftfiles[istx]=files[ifile] if (shiftfiles[istx] eq '') then begin print,' ##### sdo_stx_combine abort: no hmicont file in dir'+instx[istx] return endif endfor ; cut shiftfiles to short duration to find shifts after reform shortfiles='/tmp/shortfile_'+trim(string(indgen(nstx)))+'.fits' itmid=fix(mean(nt_arr)/2.) timearr_short=timearr_out[itmid-2:itmid+2] for istx=0,nstx-1 do begin timearr_in=starttai_arr[istx]+indgen(nt_arr[istx])*cadence_arr[istx] reformcubefile,shiftfiles[istx],shortfiles[istx],verbose=0,$ congridfactor=1./[pxratio_arr[istx],pxratio_arr[istx]],$ cutcentralx=nxcut_out,cutcentraly=nycut_out,$ smear=3,$ timearr_in=timearr_in,timearr_out=timearr_short endfor ; take mean over this short duration shortmean=fltarr(nxcut_out,nycut_out,nstx) for istx=0,nstx-1 do shortmean[*,*,istx]=getmeanfitscube(shortfiles[istx]) ; find the shifts from first images, pairwise dirs to first dir file stxshifts=fltarr(nstx,2) ; now istx=0 has shift [0,0] ; rotate first image to standard im_sdo0=shortmean[*,*,0] reformimage,im_sdo0,im_sdo0_rot,rotate=-delangle_arr[0] ; measure the other against this one after rotation to standard for istx=1,nstx-1 do begin im_sdo2=shortmean[*,*,istx] reformimage,im_sdo2,im_sdo2_rot,rotate=-delangle_arr[istx] ; find shift to put rotated 2 on rotated 0 shift_21=findimshift_rr(im_sdo0_rot,im_sdo2_rot,/subpix,$ filter=20,cropborders=1,croptriangles=1) print,' ===== '+shortfiles[istx]+' has shift to '+shortfiles[0]+' = '+$ trimd(shift_21) stxshifts[istx,0]=shift_21[0] stxshifts[istx,1]=shift_21[1] ;; ; check ;; reformimage,im_sdo2_rot,im_shifted,shift=shift_21 ;; showex,im_sdo0_rot,im_shifted ;; STOP endfor ; ----- selectfiles specified: handle only these files if (selectfiles[0] ne '') then begin files=selectfiles nfiles=n_elements(files) for ifile=0,nfiles-1 do begin ; find istx for this file istx=-1 filepath=file_dirname(files[ifile]) for iistx=0,nstx-1 do $ if (strmatch(filepath,stxnames[iistx]+'2rotsdo')) then istx=iistx if (istx eq -1) then begin print,' ##### sdo_stx_combine abort: no dirname for selectfiles file path' return endif fileonly=file_basename(files[ifile]) inselstring=cgrootname(fileonly) outselstring='_'+stxnames[istx]+'_comb' timearr_in=starttai_arr[istx]+indgen(nt_arr[istx])*cadence_arr[istx] doallfiles,stxnames[istx]+'2rotsdo',outdir,inselstring,outselstring,$ 'reformcubefile',$ congridfactor=1./[pxratio_arr[istx],pxratio_arr[istx]],$ shift=[stxshifts[istx,0],stxshifts[istx,1]],$ rotate=-delangle_arr[istx],$ cutcentralx=nxcut_out,cutcentraly=nycut_out,$ timearr_in=timearr_in,timearr_out=timearr_out,$ missingvalue=-1,splinip=splinip endfor ; end loop over selectfiles files endif ; end test sellib present ; ----- if no selectfiles then process all files if (selectfiles[0] eq '') then begin for istx=0,nstx-1 do begin ; define regular timearr_in timearr_in=starttai_arr[istx]+indgen(nt_arr[istx])*cadence_arr[istx] ; in first dir do all SDO files if (istx eq 0) then begin ; AIA files inselstring='aia*fits' outselstring='_'+stxnames[istx]+'_comb' doallfiles,stxnames[istx]+'2rotsdo',outdir,inselstring,outselstring,$ 'reformcubefile',$ congridfactor=1./[pxratio_arr[istx],pxratio_arr[istx]],$ shift=[stxshifts[istx,0],stxshifts[istx,1]],$ rotate=-delangle_arr[istx],$ cutcentralx=nxcut_out,cutcentraly=nycut_out,$ timearr_in=timearr_in,timearr_out=timearr_out,$ missingvalue=-1,splinip=splinip ; HMI files inselstring='hmi*fits' outselstring='_'+stxnames[istx]+'_comb' doallfiles,stxnames[istx]+'2rotsdo',outdir,inselstring,outselstring,$ 'reformcubefile',$ congridfactor=1./[pxratio_arr[istx],pxratio_arr[istx]],$ shift=[stxshifts[istx,0],stxshifts[istx,1]],$ rotate=-delangle_arr[istx],$ cutcentralx=nxcut_out,cutcentraly=nycut_out,$ timearr_in=timearr_in,timearr_out=timearr_out,$ missingvalue=-1,splinip=splinip endif ; end writing SDO files from first dir (?? oops re shift if smallest?) ; now in all 2rotsdo input dirs do only the non-AIA and non HMI-mag files outselstring='_comb' files=findfile(stxnames[istx]+'2rotsdo/*rot*.*') nfiles=n_elements(files) for ifile=0,nfiles-1 do begin fileonly=file_basename(files[ifile]) if (not(strmatch(fileonly,'*aia*') or strmatch(fileonly,'*hmimag*'))) $ then begin inselstring=cgrootname(fileonly) outselstring='_'+stxnames[istx]+'_comb' doallfiles,stxnames[istx]+'2rotsdo',outdir,inselstring,outselstring,$ 'reformcubefile',$ congridfactor=1./[pxratio_arr[istx],pxratio_arr[istx]],$ shift=[stxshifts[istx,0],stxshifts[istx,1]],$ rotate=-delangle_arr[istx],$ cutcentralx=nxcut_out,cutcentraly=nycut_out,$ timearr_in=timearr_in,timearr_out=timearr_out,$ missingvalue=-1,splinip=splinip endif endfor ; end loop over non-SDO files in this STX dir endfor ; end loop over all STX2rotsdo dirs endif ; end of no-selectfiles handling all files ; 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 ====================== ; Vasco SST CRISP + CROMIS cd,'/media/rutten/SSTNEW/alldata/SST/2018-06-10-vasco' stxnames=['cah','ha','cair'] outdir='sst2rotsdo' splinip=1 selectfiles=['cah2rotsdo/hmicont_rot_ip.fits',$ 'ha2rotsdo/hmicont_rot_ip.fits',$ 'cair2rotsdo/hmicont_rot_ip.fits'] ;;,$ ;; 'cah2rotsdo/cont_cah_rot_ip.fits',$ ;; 'ha2rotsdo/wb.6563.12:21:09.corrected_rot_ip.icube',$ ;; 'cair2rotsdo/wb.8542.12:21:09.corrected_rot_ip.icube'] selectfiles='' ; undo the above sdo_stx_combine,stxnames,selectfiles=selectfiles,outdir,splinip=splinip showex,$ /allsdo,sdodir='sst2rotsdo',$ /allmwf,mwfdirs='sst2rotsdo',$ mwfwavsfiles=['chromis/spectfile.3968.idlsave',$ 'crisp/spectfile.6302.idlsave',$ 'crisp/spectfile.6563.idlsave',$ 'crisp/spectfile.8542.idlsave'],$ /addlabels,/plotscatter ;; showex,'sst2rotsdo/hmicont_rot_ip_cah_comb.fits',$ ;; 'sst2rotsdo/hmicont_rot_ip_ha_comb.fits',$ ;; 'sst2rotsdo/hmicont_rot_ip_cair_comb.fits' ;; ,$ ;; 'sst2rotsdo/cont_cah_rot_ip_cah_comb.fits',$ ;; 'sst2rotsdo/wb.6563.12:21:09.corrected_rot_ip_ha_comb.icube',$ ;; 'sst2rotsdo/wb.8542.12:21:09.corrected_rot_ip_cair_comb.icube' end