; file: sdo_makeslavefitscube.pro ; init: Dec 30 2013 Rob Rutten Deil ; last: Nov 15 2020 Rob Rutten Deil ; note: change "slave" into politically correct softonom? ;+ pro sdo_makeslavefitscube,datadir,wav,refwav=refwav,$ notrack=notrack,diffdesolrot=diffdesolrot,border=border,$ nocrossalign=nocrossalign,driftsdir=driftsdir,$ clipmin=clipmin,clipmax=clipmax,no_intscale=no_intscale,$ euvanchor=euvanchor,itend=itend,verbose=verbose ; PURPOSE: ; convert level 1.5 SDO images at slave wavelength into fitsfile cube, ; with fractional-pixel derotation, cadence interpolation to the reference, ; intensity rescaling, optional crossaligning, optional derotation ; ; CALL EXAMPLE: ; see end of program file ; ; INPUTS: ; datadir = string path/dir containing dir /level2 with SDO files ; wav: string specifying level2 diagnostic, e.g. '1700', 'mag' ; ; OPTIONAL KEYWORD INPUTS: ; refwav: string specifying reference wavelength, default '171' ; notrack = 1/0: JSOC did not track rotation in im_patch cutting out ; diffdesolrot = 1/0: differentially derotate target ; border: only here to increase if necessary ; nocrossalign = 1/0: default 0 = cross-align per driftsplines ; driftsdir: path to driftsplines, default datadir+'/../driftscenter' ; clipmin: min intensity clipping in sdo_intscale.pro ; clipmax: max intensity clipping in sdo_intscale.pro ; no_intscale = 1/0: do not intscale target AIA and HMIcont intensities ; euvanchor: anchor EUV in reversed order, default ['304','mag'] ; itend: earlier stop for quicker tests (default 0 = do all) ; verbose = 1/0: print and plot shifts per timestep ; ; OUTPUTS: ; cube fits file ; ; METHOD: ; use assoc to permit larger cubes than memory size ; ; HISTORY: ; Jan 9 2014 RR: start ; Sep 12 2015 RR: use headfits_rr to accommodate Mac dir with space ; Feb 25 2016 RR: crossalign: also EUV > HMI alignment ; Jan 10 2018 RR: option no_intscale ; Jul 5 2020 RR: diffdesolrot for target, auto-increase border ;- ; answer no-parameter query if (n_params() lt 2) then begin sp,sdo_makeslavefitscube return endif ; defaults for keywords if (n_elements(refwav) eq 0) then refwav='171' if (n_elements(notrack) eq 0) then notrack=0 if (n_elements(diffdesolrot) eq 0) then diffdesolrot=0 if (n_elements(border) eq 0) then border=15 if (n_elements(nocrossalign) eq 0) then nocrossalign=0 if (n_elements(driftsdir) eq 0) then driftsdir=datadir+'/../driftscenter' if (n_elements(clipmin) eq 0) then clipmin=0 if (n_elements(clipmax) eq 0) then clipmax=0 if (n_elements(no_intscale) eq 0) then no_intscale=0 if (n_elements(euvanchor) eq 0) then euvanchor=['304','mag'] if (n_elements(itend) eq 0) then itend=0 if (n_elements(verbose) eq 0) then verbose=0 ; define datadir subdirs level2dir=datadir+'/level2' cubesdir=datadir+'/cubes' ; make cubes dirs and plots subdir if not yet existing spawn,'mkdir -p '+cubesdir if (verbose ne 0) then spawn,'mkdir -p '+cubesdir+'/plots_align/' ; define permitted level2 SDO (AIA + HMI) filename strings wavs=['94','131','171','193','211','304','335','1600','1700',$ ; AIA 'mag','cont','dop'] ; HMI nwavs=n_elements(wavs) fitscubenames=['aia94.fits','aia131.fits','aia171.fits',$ 'aia193.fits','aia211.fits','aia304.fits',$ 'aia335.fits','aia1600.fits','aia1700.fits',$ 'hmimag.fits','hmicont.fits','hmidop.fits'] ; set channel names for plots channels=wavs for i=0,8 do channels[i]='AIA '+wavs[i] channels[9]='HMI magnetogram' channels[10]='HMI continuum' channels[11]='HMI Dopplergram' ; check wav specification and get iwav iwav=-1 for iw=0,nwavs-1 do if (wavs[iw] eq wav) then iwav=iw if (iwav eq -1 ) then begin print,' ===== ERROR: wav invalid = ',wav print,' ===== valid: ',wavs return endif ; get iref = index of reference (and check it is there) iref=-1 for iw=0,nwavs-1 do if (wavs[iw] eq refwav) then iref=iw if (iref eq -1 ) then begin print,' ===== ERROR: refwav invalid = ',wav print,' ===== valid: ',wavs return endif refcube=fitscubenames[iref] ; uncompress level2 files if necessary print,' ----- uncompress level2 wav files if needed (may take long)' spawn,'funpack -D '+datadir+'/level2/*'+wav+'.fits.fz' ; get file list files=findfile(datadir+'/level2/*'+wav+'.fits') nfiles=n_elements(files) ; get refcube size (was cropped by border parameter) refheader=headfits_rr(cubesdir+'/'+refcube,/silent) nxref=fxpar(refheader,'naxis1') nyref=fxpar(refheader,'naxis2') cenxpxref=nxref/2.-0.5 ; center of FOV center px in IDL no-finger count cenypxref=nyref/2.-0.5 ntref=fxpar(refheader,'naxis3') if (itend eq 0) then itend=ntref-1 ; get reference timing array starttim=fxpar(refheader,'starttim') tairef=anytim2tai(starttim) cadence=fxpar(refheader,'cadence') timarr=indgen(ntref)*cadence+tairef ; get reference FOV location data cdelt1ref=fxpar(refheader,'cdelt1') cdelt2ref=fxpar(refheader,'cdelt2') xcenarcref=fxpar(refheader,'xcen') ; FOV center solar X [arcs from sun center] ycenarcref=fxpar(refheader,'ycen') crota2ref=fxpar(refheader,'crota2') rsun_obs=fxpar(refheader,'rsun_obs') ; define outheader with timing and location parameters for future use ; set basic parameters bigendian=1 ;RR my Toshiba laptop bitpix=16l ;RR also for HMI cubes mkhdr,outheader,abs(bitpix)/8,[nxref,nyref,itend+1] sxaddpar,outheader,'starttim',starttim,'first image mid-exposure' sxaddpar,outheader,'cadence',cadence,'image cadence [s]' sxaddpar,outheader,'cdelt1',cdelt1ref,'x-scale [arcsec/px]' sxaddpar,outheader,'cdelt2',cdelt2ref,'y-scale [arcsec/px]' sxaddpar,outheader,'crota2',crota2ref,$ 'rotation angle degrees NESW = CCW' sxaddpar,outheader,'rsun_obs',rsun_obs,'solar radius [arcsec]' sxaddpar,outheader,'xcen',xcenarcref,$ 'FOV center X [arcsec from sun center]' sxaddpar,outheader,'ycen',ycenarcref,$ 'FOV center Y [arcsec from sun center]' sxaddpar,outheader,'channel',channels[iwav],$ 'SDO diagnostic name [string]' ; get size outheader sizeoutheader=size(outheader) outheadersize=(1+fix(sizeoutheader[1]/36.))*2880 ;RR fits header = Nx36 "card images" = Nx2880 bytes ; open output cube file for assoc, write header get_lun, unit_out if (bigendian) then openw,unit_out,cubesdir+'/'+fitscubenames[iwav],$ /swap_if_little_endian $ else openw,unit_out,cubesdir+'/'+fitscubenames[iwav] outassoc=assoc(unit_out,intarr(nxref,nyref),outheadersize) rec=assoc(unit_out, bytarr(outheadersize)) rec[0]=byte(outheader) ; read first image as initial image0 shiftx=fltarr(nfiles) shifty=fltarr(nfiles) read_sdo,files[0],index0,image0,/uncomp_delete sizeim0=size(image0) nx=sizeim0[1] ny=sizeim0[2] ; =============== first image ; rescale first image image0=sdo_intscale(index0,image0,wav,$ clipmin=clipmin,clipmax=clipmax,no_intscale=no_intscale) time0=anytim2tai(index0.t_obs) missing=avg(image0) ; if crossalign get AIA shift against HMI for first image sdoshift=[0,0] if (iwav lt 9 and nocrossalign ne 1) then $ sdoshift=sdo_getsumsplineshift(wav,time0,driftsdir,euvanchor=euvanchor) if (verbose ne 0) then $ print,' ======== wav, sdoshift(it=0): ',wav,trim(sdoshift) ; reposition first image sdo_shiftjsoc,image0,index0,$ tairef,xcenarcref,ycenarcref,cenxpxref,cenypxref,sdoshift,$ outimage,outshift,notrack=notrack,diffdesolrot=diffdesolrot image0=outimage shiftx[0]=outshift[0] shifty[0]=outshift[1] ; ============= second image ; read second image as initial image1 read_sdo,files[1],index1,image1,/uncomp_delete time1=anytim2tai(index1.t_obs) ; rescale second image image1=sdo_intscale(index1,image1,wav,$ clipmin=clipmin,clipmax=clipmax,no_intscale=no_intscale) ; if crossalign get AIA shift against HMI for second image if (iwav lt 9 and nocrossalign ne 1) then $ sdoshift=sdo_getsumsplineshift(wav,time1,driftsdir) ; reposition second image sdo_shiftjsoc,image1,index1,$ tairef,xcenarcref,ycenarcref,cenxpxref,cenypxref,sdoshift,$ outimage,outshift,notrack=notrack,diffdesolrot=diffdesolrot image1=outimage shiftx[1]=outshift[0] shifty[1]=outshift[1] ; start loop over timesteps itfile=1 for it=0,itend do begin ; ============== next image ; get a new image if image1 is still before the new it time if (time1 le timarr[it] and itfile lt nfiles-1) then begin ; previous image1 now becomes image0 image0=image1 time0=time1 ; read new image1 itfile=itfile+1 read_sdo,files[itfile],index1,image1,/uncomp_delete time1=anytim2tai(index1.t_obs) ; intscale new image image1=sdo_intscale(index1,image1,wav,clipmin=clipmin,clipmax=clipmax,$ no_intscale=no_intscale) ; if crossalign then get AIA shift against HMI for this new image if (iwav lt 9 and nocrossalign ne 1) then $ sdoshift=sdo_getsumsplineshift(wav,time1,driftsdir) ; reposition new image sdo_shiftjsoc,image1,index1,$ tairef,xcenarcref,ycenarcref,cenxpxref,cenypxref,sdoshift,$ outimage,outshift,notrack=notrack,diffdesolrot=diffdesolrot image1=outimage shiftx[itfile]=outshift[0] shifty[itfile]=outshift[1] endif ; end getting next image for straddle pair ; interpolate current straddle pair to refwav timing time step ;RR pairwise rotation and drift corrected so only evolutionary differences p=[[[image0]],[[image1]]] timefrac=(timarr[it]-time0)/(time1-time0) imageit=interpolate(p,timefrac) ;RR bilinear, cubic no difference for 2) ; print info if (verbose ne 0) then begin if (it eq 0) then $ print,' it, itfile, time, dtairef, dt1, tfrac, shifts x,y' print,' ',ntostr(it),' ',ntostr(itfile),' ',$ ntostr([timarr[it]-tairef,time0-tairef,time1-tairef,timefrac,$ shiftx[itfile],shifty[itfile]],format='(F12.2)') endif ; check border and auto-increase its value when too small if (nx lt nxref or ny lt nyref) then begin border=border+15 print,' ##### nx, ny = '+ntostr(nx)+' '+ntostr(ny)+$ ' but nxref, nyref = '+ntostr(nxref)+' '+ntostr(nyref) print,' if you type ".cont" sdo_maketargetfitscubes will run with ' print,' border ='+trimd(border)+' but no mpegs or targetmovie' STOP spawn,' cd '+datadir sdo_maketargetfitscubes,targetdir=targetdir,refwav=refwav,$ notrack=notrack,diffdesolrot=diffdesolrot,border=border,$ nocrossalign=nocrossalign,driftsdir=driftsdir,$ clipmin=clipmin,clipmax=clipmax,no_intscale=no_intscale RETALL endif ; write interpolated image cropped to reference size cutim=imageit[0:nxref-1,0:nyref-1] ; cuts over lowerleft borders outassoc[it]=cutim ; check ;; if (it eq 0) then print,' ##### nx,ny,nxref,nyref = ',nx,ny,nxref,nyref endfor ; end loop over it time steps ; free the output file free_lun,unit_out ; plot the image-center shifts if (verbose ne 0) then begin shiftx=shiftx[0:itfile] shifty=shifty[0:itfile] psfilename=cubesdir+'/plots_align/shifts-'+fitscubenames[iwav] psfilename=str_replace(psfilename,'.fits','.eps') openpsplot,psfilename,thick=2,fontsize=9,xsize=8.8,ysize=8.8 yrange=[min([min(shiftx),min(shifty)]),$ max([max(shiftx),max(shifty)])] yrange=[yrange[0]-0.1*(yrange[1]-yrange[0]),$ yrange[1]+0.1*(yrange[1]-yrange[0])] plot,indgen(itfile),shiftx[0:itfile],$ xrange=[0,itfile],yrange=yrange,xstyle=1,ystyle=1,$ xtitle='time step',ytitle='shift in px',$ position=[0.2,0.1,0.95,0.95] oplot,indgen(itfile),shifty[0:itfile],linestyle=2 itlabel=(itfile)/2 labeloffset=(yrange[1]-yrange[0])*0.02 xyouts,itlabel,shiftx[itlabel]+labeloffset,'x' xyouts,itlabel,shifty[itlabel]+labeloffset,'y' xyouts,itfile/8,yrange[1]-0.1*(yrange[1]-yrange[0]),'shifts '+fitscubenames[iwav],$ charsize=1.2 closepsplot,psfilename,opengv=verbose ; print success print,' === wrote ',cubesdir+'/'+fitscubenames[iwav] endif ; end of shifts plotting end ; =============== main for testing per IDLWAVE H-c ======================== ;; datadir='/home/rutten/data/SDO/2014-06-14-small/center' ;; datadir='/home/rutten/data/SDO/2014-06-14-small/target' datadir='/home/rutten/data/SDO/2019-11-11-transit/midpointfull/target' wav='cont' ;; diffdesolrot=0 ; center diffdesolrot=1 ; target ;; nocrossalign=1 ; center nocrossalign=0 ; target sdo_makeslavefitscube,datadir,wav,$ diffdesolrot=diffdesolrot,$ nocrossalign=nocrossalign,driftsdir=driftsdir,/verbose ; check showex,datadir+'/cubes/hmi'+wav+'.fits' ; ########## adapt > hmi or aia end